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Abstract 

We consider the 2+1 and 3+1 scalar wave equations reduced via a helical Killing 
field, respectively referred to as the 2-dimensional and 3-dimensional helically re- 
duced wave equation (HRWE). The HRWE serves as the fundamental model for 
the mixed-type PDE arising in the periodic standing wave (PSW) approximation 
to binary inspiral. We present a method for solving the equation based on domain 
decomposition and spectral approximation. Beyond describing such a numerical 
method for solving strictly linear HRWE, we also present results for a nonlinear 
scalar model of binary inspiral. The PSW approximation has already been theoreti- 
cally and numerically studied in the context of the post-Minkowskian gravitational 
field, with numerical simulations carried out via the "eigenspectral method." Despite 
its name, the eigenspectral technique does feature a finite-difference component, and 
is lower-order accurate. We intend to apply the numerical method described here 
to the theoretically well-developed post-Minkowski PSW formalism with the twin 
goals of spectral accuracy and the coordinate flexibility afforded by global spectral 
interpolation. 
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1 Introduction and preliminaries 



1.1 Periodic standing wave (PSW) approximation 

The development of gravitational wave detectors has spurred interest in the 
binary inspiral of compact astrophysical objects, in particular black holes. 
The challenge of computationally solving Einstein's equations for such a sys- 
tem has become the focus of many groups throughout the world both for 
its importance in gravitational wave astronomy and for its role in advanc- 
ing the understanding of highly dynamical strongly curved spacetimes. Very 
recentlyflBSSSi,^ a number of techniques have been found to stabilize 
codes, so that computational evolution of Einstein's equations can track the 
last few orbits of binary inspiral. 

There is, of course, more to the problem than the last few orbits. For orbit- 
ing objects with mass M and separation a, the characteristic measure of the 
nonlinearity of their gravitational interaction is GMja(? where G^c are the 
gravitational constant and the speed of light. When this measure is small, 
the dynamics and the generation of gravitational waves can be found using 
post-Newtonian methods, an analytic approach in which Einstein's theory is, 
in effect, expanded treating GMja(? as a perturbation parameter. 

With methods available for the small-GM/ac^ early stage, and the large- 
GMj a(? last few orbits, what remains needed is an effective method for treat- 
ing the intermediate epoch, the stage of inspiral in which GMj a(? is too large 
for post-Newtonian methods, but in which many orbits remain. When more 
than a few orbits still remain, an accurate numerical evolution of Einstein's 
equations will be too computationally expensive, at least in the near future. 

The periodic standing wave (PSW) approach is an approximation scheme for 
dealing effectively with this intermediate inspiral epoch. The scheme is based 
on the fact that during this epoch the orbiting inspiral is very much orbit, and 
only slightly inspiral. That is, the change in the orbital radius is small for each 
orbit. (This is, in fact, a criterion for many orbits to remain, and for accurate 
computational evolution to be daunting.) In the PSW approach the motion 
of the sources and the fields are assumed to be helically symmetric, that is, 
all quantities are rigidly rotating in the sense that a change in time by At is 
equivalent to a change in azimuthal angle according to A(y9 — > — f2At, where 
r2 is a constant, the angular velocity with which the fields rigidly rotate. 

The imposition of this helical symmetry vastly changes the nature of the math- 
ematical and computational problem. Prior to helical reduction, the problem is 
that of evolving forward in time a hyperbolic problem (more precisely, a prob- 
lem that can be cast in hyperbolic form). The imposition of helical symmetry 
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reduces by one the number of independent variables and, more important, 
changes the problem from a hyperbolic problem to a mixed one, a problem 
with a region of the manifold (near the rotation axis) in which the equations 
are elliptic, and an outer region in which the equations are hyperbolic. The 
boundary conditions for this problem are also unusual. One can have the pres- 
ence of a source represented by inner Dirichlet boundary conditions on two 
small topological spheres just outside the location of the sources, the compact 
astrophysical objects. Alternatively, one can include the objects as explicit in- 
homogeneities in the equations. The other boundary conditions on the problem 
represent the radiation in the distant wave zone, and require some discussion. 



An important feature of this problem in general relativity is the conservation 
of the total energy of the system. If energy is leaving in the form of outgoing 
gravitational radiation, then the orbital motion cannot be helically symmetric; 
the radius must decrease. For fields other than gravity one could invoke a 
force, some deux ex machina, to keep the orbits unchanging, and have that 
force not couple to the field being studied. This certainly can be, and has been 
done in model problemsfs], 0, [lO, 11, 12] but, in principle, cannot be done in 



general relativity. In Einstein's gravitation all forces couple to gravity. This is 
dealt with by computing a standing wave solution of the helically symmetric 
problem, that is, by imposing standing wave outer boundary conditions. From 
the standing wave exact solution, an approximate solution is then extracted 
to the physical problem of outgoing waves. 



By "standing waves" in a linear theory, we mean here the average of the 
outgoing wave solution and the ingoing wave solution. No clear meaning exists 
for standing waves in a nonlinear theory. To define our standing waves, we 
choose an iterative solution technique for the nonlinear problem in which an 
average of ingoing and outgoing solutions is taken at each step of iterationll] 
We emphasize here that in each iteration of the PSW problem it is an outgoing 
(or, trivially different, ingoing) problem that is solved. The present paper will 
therefore focus on the details of the efficient computation of a problem with 
outgoing radiation boundary conditions. 

Hyperbolic problems with Dirichlet boundary conditions are known not to be 
well posed. We take a pragmatic approach to whether our mixed-type problems 
with radiative conditions are well posed. For one thing, the problem arises in 
what would appear to be a physically well-specified problem; heuristically 
the problems of nonuniqueness for the Dirichlet problem are removed by the 
radiation boundary conditions. For another thing, no fundamental instability 
has been encountered in seeking numerical solutions of our problems. Further 



Another extension of "standing wave" to nonlinear theories has been discussed in 
the Uterature[^, using the minimum amplitude of the multipole in every multipole 
mode. 
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indirect evidence that the problem is well posed can be found in the work 
of Torre, who has shown that a closely related mixed- type linear problem is 
well-posed if the boundary conditions are of an admissible type that includes 
Sommerfeld outgoing conditions 13|, [14 . 



The strategy of the PSW approach is to solve the helically symmetric binary 
problem computationally, but otherwise without approximation. From that 
"exact" helically symmetric solution, an approximation is extracted for the 
physical problem. The extraction of an outgoing approximation is tantamount 
to treating the nonlinear standing wave solution as if it were an average of the 
ingoing and outgoing solutions. (For details of extraction, see lij.) The rea- 
son that this is an excellent approximation (as is confirmed by computations 
with model problems) is that the regions of the physical manifold in which 
nonlinearities are strong are those near the sources, and in these regions the 
solution is extremely insensitive to the boundary conditions (ingoing, outgo- 
ing, standing wave) in the distant weak field wave region. Since the ingoing 
and outgoing solutions are nearly identical in this region they can be averaged 
(or simply replaced by either the ingoing or the outgoing solution). Where 
the ingoing/outgoing solutions are very different, the theory is approximately 
linear and hence again the ingoing and outgoing fields can be averaged. This 
feature, the separation of the region of nonlinearity and the region of waves, 
is closely related to the argument, given below, that a multidomain spectral 
method should have advantages for PSW type problems even beyond the ad- 
vantages it has demonstrated for purely elliptic initial-value problems in the 



work of Pfeiffer et a l. |l5l . Il6l |. or its use in evolution codes by the Caltech- 
Cornell collaboration 17 . 



Successful solutions of the PSW approximation will serve a variety of purposes. 
As discussed above, it will provide the bridge between the post-Newtonian 
methods and numerical evolution for binary orbits; it will provide near optimal 
starting points for the numerical evolution of the last few orbits; it may give a 
useful testbed for studying radiation reaction; solution of the PSW standing 
wave problem will give a new class of solutions to Einstein's theory. 



Work on PSW computations has already taken several steps using the "eigen- 
spectral method," an approximate numerical method based on coordinates 
well adapted to the geometry both close to the sources and in the radiation 
region. This method, in effect, keeps only the features of the solution that are 
most important to the structure of the near-source fields and to the radiation, 
and has been used so far to study nonlinear toy models|8l,0,[i3, linearized gen- 
eral relativity 18|] , and the post-Minkowski extension of linearized gravity 19 . 
This method is remarkable for its simplicity in giving approximate solutions, 
but is ill suited to the high accuracy needed for several purposes. In addition 
the eigenspectral technique uses a finite-difference component, and provides 
numerical answers on a specific grid of points in the space of the physical prob- 



4 



lem. Its results, therefore, require a difficult interpolation from a nonuniform 
grid if they are to be compared with other results or used as initial conditions 
for evolution codes. The multidomain spectral method we will describe has 
neither shortcoming. 

To minimize the complexity of issues not directly relevant to the multidomain 
spectral method, we choose as a specific target problem, more-or-less the prob- 
lem of Refs. js], M, [13] , that of a nonlinear scalar field. The physical problem, 
of course, is one with three spatial coordinates, but for simplicity both of ex- 
position and of computation, we work here with the two-dimensional helically 
reduced wave equation, hereafter 2d HRWE. The Cauchy problem for this 
physical model would involve three independent variables (two spatial, one 
time), but the helical reduction means we are solving on a manifold of two 
dimensions. We will also, in Sec. 4, present some preliminary results for the 3d 
HRWE, and in the concluding section will discuss application of the method 
to the 3d HRWE problem. The physical picture is that of two compact scalar 
charges moving on circular orbits of radius a, with angular velocity fl, and 
emitting helically symmetric outgoing waves. The wave speed c will be taken 
to be unity. The scalar field satisfies a 2d HRWE of the form 



on the region outside the two sources and inside a large radius circle (spherical 
surface for 3d HRWE) on which outgoing conditions are imposed. Here L is 
a mixed-type linear HRWE operator (the helically reduced d'Alembertian), h 
is a nonlinear function of ip, and the i] parameter controls the strength of the 
nonlinearity. The inhomogeneity g can be used to represent explicit sources. 
In our 2d HRWE development below we will use inner boundary conditions, 
not explicit sources, so g will be zero; a nonzero g will be used only in the 
consideration of the 3d HRWE in the concluding section. 

We make one more simplifying choice. We take the solution to have the sym- 
metry ipi—x, —y) = —ipix, y), a choice equivalent to taking the sources to have 
opposite scalar charge. By setting the total charge equal to zero, we eliminate 
the bothersome logarithmic monopole that would diverge at large distances. 
Of course, for gravitation theory there is only charge of a single sign and the 
monopole cannot be set to zero, but the real physical problem is that for 
three spatial dimensions, in which the monopole falls off with distance from 
the sources, and is an acceptable complication. 

1.2 Multidomain spectral method 

The "two center" domain of the 2d HRWE problem is illustrated in Fig. 1. 
This figure shows how two annular domains, H and A along with eight rect- 



Lip + r]h{ip) = g , 
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Fig. 1. The 2d HRWE problem and its domain decomposition. The circles 
that are the outer boundaries of the darkened source regions carry Dirichlet bound- 
ary conditions representing the imprint of the sources. Region H, between the two 
smah circles concentric with the source, is the inner annular domain. Domain A 
between two large solid circles is the outer annular domain. The dashed circle in 
domain A is the light circle. 

angular domains, fully cover the physical problem. Note that the symme- 
try ip{—x, —y) = —ip{x,y) means that annulus H represents the information 
around both the source on the right, and that on the left. Similarly, rectan- 
gular domains 1,2,3,7,8 carry the information about the solution in symmetry 
related regions on the left side of the physical problem. Using the aforemen- 
tioned symmetry, we could do away with region 6, which is equivalent to 
region 4. Nevertheless, we have chosen to keep region 6 so that our code can 
be tested on elliptic problems which, when posed on the inner region spanned 
by the rectangles, need not be of definite parity. The annular domain A ex- 
tends between the two solid circles in Fig. 1. At its outer boundary we impose 
Sommerfeld-like outgoing radiation boundary conditions to be described be- 
low. For most radiation conditions to be applicable, the outer boundary of 
A must be at least several wavelengths away from the sources. For a typical 
choice of Q this means that the radius of the outer boundary must be 10 times 
or more larger than the separation between the source regions. The dashed 
circle in Fig. 1 is the "light circle," the boundary between the inner region in 
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3,7 



2,8 







Fig. 2. Domain decomposition for the 3d HRWE Problem. Here the inner 
spherical shell and outer spherical shell, respectively corresponding to the annuli H 
and A in Fig. 1, are not shown. The remaining domains (all lying within the elliptic 
region) are shown in an exploded format for emphasis. 

which the operator L is elliptic, and the outer region in which L is hyperbolic. 
Typically the radius of the the light circle will be at least several times larger 
than the separation of the centers of the source domains. 



For a 3d problem, the analogous two-center domain is actually composed of 
fewer elements than in the equivalent 2d problem. This is one of the reasons 
we believe that a multidomain spectral method is well-suited for the compu- 
tations needed in PSW approximation of binary inspiral. In lieu of a detailed 
description of the domain decomposition we envision for 3d work, we offer the 
picture in Fig. 2. Note that cylindrical shell labelled 2,8, for example, corre- 
sponds to both domain 2 and domain 8 in Fig. 1. Indeed, the 2d figure could 
alternatively be viewed as a cross section of a 3d scenario. 

The ten subdomains H,l — 8, A, along with the symmetry ■?/;(— x, —y) = 
—ip{x,y), contain the complete information about ip. Below we describe a 
spectral method in which the solution for ip in each region is considered as 
an expansion in terms of appropriate basis functions. Moreover, we use the 
term spectral in the truest sense, since the unknowns we solve for in our nu- 
merical approach are in fact the coefficients of the basis functions. (In a later 
section, while considering the 3d HRWE we will also describe results based 
in part on a pseudospectral method in which the unknowns are point-values 
on a spectral grid.) A crucial feature of our numerical method is to separate 
the physical problem into subdomains so that the type changing nature of the 
HRWE is then confined within a single subdomain, the outer annulus A. The 
inner region is then purely elliptical and, despite its two center topology, this 
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inner part of the manifold is amenable to the standard spectral techniques 
associated with elliptic equations 15|, [16 . 



1.3 Brief review of mixed-type problems 



In order to place our work in context, we offer some remarks on type-changing 
PDE in general. Perhaps the most notable examples arise in the mathemat- 
ical description of transonic flow, and in particular flow over an air foil, a 
scenario for which subsonic and supersonic regions are respectively described 



by elliptic and hyperbolic equations. The Frankl-Chaplygin equation [i^ 
K{y)Uxx + Uyy = (the choice K{y) = y determines the familiar Tricomi 
equation 2^) serves as a fundamental model of this behavior, and Morawetz 
carried out early and fundamental analytic studies 2^ of it and correspond- 
ing first order systems (see also 2^). The first truly successful method for 
numerically calculating transonic flow was put forth in the seminal paper 24 
by Murman and Cole. They adopted a relaxation-like method along with 
type-dependent finite difference stencils. For an account of the impact of 
the Murman-Cole method in aerodynamics, and a review of more modern 



CFD approaches towards transonic flow, see the historical perspective of [25 



While the aerodynamic scenario is the most widely known, special examples 
of mixed-type equations arise in fields as diverse as plasma physics [1^ and 
windshield design 27|. An early work most closely related to our own is Chao- 
hao's examination 28|] of amplifying spiral wave solutions to the 2-1-1 wave 
equation. Numerical methods for mixed-type problems tend to be equation 
specific. The Murman-Cole technique, for example, would seem to have no 
application for our problem. 

Motivated by mixed-type problems, in the late 1950s Friedrichs initiated a 
program 29|] for analyzing a wide class of boundary value problems based on 



operators whose symmetric part is positive definite and which obey certain ad- 
missible boundary conditions. Such symmetric positive systems include PDE 
of hyperbolic, elliptic, and mixed type. A lucid history of the Friedrichs pro- 
gram from both theoretical and numerical perspectives is given by Jensen |30|. 
Rather early on, Katsanis developed a numerical method for solving Friedrichs 



systems [3l[. Starting with a Friedrichs system, he applied Green's theorem 
in a generic cell, and approximated the resulting integral equation. At the 
discrete level the approach faithfully mimicked both both the symmetric pos- 
itive aspect of the operators as well as admissibility of boundary conditions, 
and Katsanis went on to numerically examine the Tricomi problem 32|. Al- 
though geometrically flexible, the Katsanis method is low-order accurate. A 



33 



finite-difference method for Friedrichs system was also outlined by Liu 
Recent work towards numerically solving Friedrichs systems has drawn on 
the powerful framework of discontinuous Galerkin methods MS. Although 
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Torre has shown that the 2d HRWE on a disk can be cast into a first-order 
Friedrichs system (and, indeed, the aforementioned work by Chaohao also 
made connections with Friedrichs theory), we have nevertheless chosen a clas- 
sic multidomain spectral method over a discontinuous Galerkin method. This 
would seem appropriate given the relatively simple geometry of our problem 
and the expected smoothness of the solutions we seek. 



2 Outer annulus and outer boundary conditions 



We begin by discussing the 2d HRWE on the outer annulus, labeled A in Fig. 1; 
with appropriate changes this discussion will also apply to the 3d HRWE on an 
outer spherical shell. As the relevant PDEs are type-changing on them, these 
are the most interesting subdomains. In addition, our discussion here will 
supply some of the details of the HRWE and will exhibit its mathematical 
features, including its change of type. 

In this section, we will describe the radiative outer boundary conditions on A 
imposed at a large r = rmax = R- We shall also speak of an inner boundary 
condition for A, although when A is considered as a subdomain it is not asso- 
ciated with a true inner boundary condition. The philosophy here is that we 
must understand the relevant boundary value problem on each subdomain to 
have all subdomains successfully "glued" together. Inner boundary conditions 
on A are imposed at a radial value r = rmin = e . Note that e is not small. 
For the decomposition illustrated in Fig. 1, in fact, e must be greater than the 
orbital radius. 



2.1 2d HRWE 



The linear 2d HRWE of Whelan, Krivan, and Price starts with —dfif) + 
'W'^ip = g and, in terms of polar coordinates r, 0, the helical reduction requires 
that ip he a. function only of r and if = (p — Qt. The 2d HRWE then takes the 
form 

r dr ^ dr^ ^ r^'^^ ^ dip"^ ' V^) ' ( ) 

where g{r, if) is a (/^-periodic source term, ({r) = 1 — r^f2^, and the coordinate 
ranges are 

e <r <R, 0< If <27T. (3) 

Equation (2) is clearly elliptic for r < and hyperbolic for r > One 
boundary value problem would be to seek solutions to this equation which 
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obey the boundary conditions 



' dip ^ dip ip ' 
dr dip 2r 



(4) 



R 



here with (/^-periodic functions a {(f) and /3{(f). This boundary value problem 
for annulus A is equivalent to the punctured disk examined by Torre [13|]. His 
proof of the well-posedness of this problem assumes that the radial endpoint 
e lies in the elliptic region, as do we by requiring e < Torre's proof places 
no restriction on R. He allows for R to lie in the hyperbolic region, the elliptic 
region, or even on the light cylinder r = \^\~^- In practice of course, R is large 
and lies in the hyperbolic region, as we will assume. 

Fourier transformation of (2) yields 




n 



- ^C{r)ipn = (Jnir) 



(5) 



and for the corresponding boundary conditions on the mode ipn{f) we have 



1 ( \ - 1^'^^ • n ? , ' 
ipn[e) = an, I - milipn + 



/3n 



(6) 



R 



For now we will assume n 7^ 0, and consider the zero-mode case separately 
later on. As an alternative we may also impose exact outgoing-wave boundary 
conditions, enforced n-by-n via 



— milipn + — 

dr 2r 



R 



1 

R 



nQR 



Vn{nVtR) 



(7) 



where 



VJz) 



exp 



is set up to satisfy Vu{z) ~ 1 as 2 — >• 00. Here Hl~^\z) is the first cylindrical 
Hankel function. The "frequency-domain kernel" 



nVLR 



V:.innR) 
VninQR) 



(9) 



can be computed as a continued fraction via Steed's algorithm[35|. Similar 
kernels appear in studies of radiation boundary conditions for time-domain 
wave propagation, for example in Refs. [36|, l37|, l38|, 1391 ]. 
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2.2 3d HRWE 



For the 3d case we will mostly use the same symbols as for the 2d case. The 
3d HRWE of Andrade et al.^ is 



1 d 



'A 



52 



^2 Qj, 



dr 



(10) 



where g{r, 6, ip) is again a yp-periodic source term, As^ is the unit two-sphere 
Laplacian, and the coordinate ranges now are 



e<r < R, 0<e <n, 0<^<27r. 



'Ill 



Equation (10) is elliptic for rsin6' < \Q\ ^ and hyperbolic for rsin^ > \Q\ ^. 
We seek solutions to this equation which obey the boundary conditions 



f dtp ^dtp tp' 
\ dr dip r 



(12) 



We consider the exact outgoing-wave conditions below. As before, we will 
assume that e < \^\~^, and explain why in a moment. 



Spherical harmonic transformation of (10) yields 



1 d 

I r" 

J.2 



em 



dr 



+ 1] 



'<Pem = km{r) , 



(13) 



and we will now take 

ipim{e) = aim, 



— imilipem H 

dr r 



em 



(14) 



R 



as the corresponding boundary conditions. We will assume m 7^ 0, and con- 
sider the m = case separately later on. As an alternative we may also 
impose exact outgoing-wave boundary conditions. Enforced mode-by-mode, 
they have a form similar to Eq. (7), 



'dip. 



dr 



im^liptrii + 



■ipem 



R 



1 

R 



rriQR 



v;^y,imnR) 

Ve+i/2{mQR) 



(15) 



To cast the radial equation stemming from the 3d HRWE into a form which 
resembles the radial equation stemming from the 2d HRWE, we substitute 

'^Pem = ieml\fr in (13), thereby finding 
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Boundary condition 


P 


Q 


Sommerfeld on tpn 
Sommerfeld on -0^^ 
Exact on -0^ 
Exact on ^l>im 


ns IH 

nVtR + lm.Vn{nVtR) 
mQR + Imu£_|_i/2("i^-R) 


1 

2 
1 

^ — Ref„(nOi2) 
1 — Rev£_|_i/2("^f^-R) 



Table 1 

Outer boundary conditions for the 2d and 3d HRWE. 



where 



In terms of uJe„ 
(16) when r = 



C(r) 



(17) 



(2£+l)2- 

2mQ/{2i + 1) 7^ 0, clearly something special occurs for 

< \fl\, so that 



here with i,m ^ 0. Note that \uj£„ 
\n\-^ < I ^- For this reason, we have required e < ^ above, ensuring 
that r = e lies within each mode's individual "elliptic region." 



2.3 General form for the outer boundary conditions 



From Eq. (5) we pass to a trigonometric rather than exponential representation 
of the transform ipn{i^) by writing 



V;*(r) exp(-inv9) + ipnir) exp{mip) (18) 
= 4^n{r) + V^*(r) cos{nip) + i ?/'„(r) - ?/'*(r) sin{nip) 
= Un{r) cosinip) + Wnir) sin(n<y9), 

with a similar splitting possible for tpimi'f'), or C,em.{f') if preferred. 

At r = R all boundary conditions, whether exact or some incarnation of 
Sommerfeld, may be expressed as follows: 

Rw'^iR) + pu^{R) + qWn{R) = 0, Ru'^iR) - pw^iR) + gu„(i?) = 0. (19) 

We list the possibilities considered so far in Table 1. Other choices of the form 
(19) are of course possible. For the exact conditions listed in Table 1, we have 
made use of the notation 

Vliz) 



Vy{z) = Z 



VJz) 



(20) 



with the obvious notation for real and imaginary parts. 



Turning now to the rather more delicate zero-modes (n = for 2d and m = 
for 3d), we note that for 2; — > the Hankel function H^^'^^z) is proportional 



12 



to z " for z/ 7^ and to log z for z/ = 0, so that in either case 

zHi-^^\z)/Hl'^\z)^-u, (21) 

and hence from the the definition in Eq. (8) we have 

v,{0) = ^-iy. (22) 

From this our key result follows. For the n = homogeneous case, Eq. (5) has 
solutions ipo{r) = c or ipo{r) oc logr. With Vo{0) = |, we find p = = g for the 
exact outgoing conditions as stated in (19). So this boundary condition rules 
out the logr solution. For the m = homogeneous case, Eq. (13) has solutions 
tjj£o OC and ^jJ£o oc r~^^~^^\ Now f £+1/2(0) = —i, and p = 0, g = £ + 1 for the 
exact outgoing boundary condition. In this case the boundary condition rules 
out the solution. 



3 Sparse spectral approximation of the 2d HRWE 



This section describes how we use spectral methods to numerically solve the 
2d HRWE on the two center domain shown in Fig. 1. Our method will ex- 
ploit the spectral integration preconditioning (IPC) proposed by Coutsias, 



Hagstrom, Hesthaven, and Torres [40l|. a general-purpose spectral method for 
solving ODEs and PDEs. Developing IPC for general orthogonal polynomi- 
als, they mostly considered it in the context of ODEs. In this context they 
presented a detailed theoretical analysis of both conditioning and convergence 
(see also [iH). Although they also carefully outlined how to apply the method 
in higher dimensions (with several illuminating two-dimensional examples), 
they did not explore conditioning issues for PDEs. 

We shall follow the IPC approach for the following reasons. Foremost, it is 
a direct recipe for applying spectral methods in a PDE setting, affords a 
straightforward way of treating both boundary conditions and the "gluing" 
of subdomains. Moreover, IPC allows us simultaneously to achieve a sparse 
banded matrix representation of the 2d HRWE on all subdomains. While 
we question whether a fully 3d multidomain PSW problem can be treated 
exclusively with spectral methods based on IPC (further comments to follow), 
we believe the method is well suited for handling the 2d HRWE on the outer 
annulus A (or the 3d HRWE on an outer spherical shell). The mixed-type 
boundary value problem on this subdomain introduces several novel features, 
and much of our analysis below centers on associated conditioning issues. In 
short, the IPC approach offers us a quick and easy way to test the fundamental 
idea of using a multidomain spectral method to solve the 2d HRWE, and has 
promise for at least one key aspect of true 3d PSW problems. 
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In following the IPC approach we are able to explicitly form the matrix which 
represents the 2d HRWE on the two center domain, and we subsequently use 
Gaussian elimination (as embodied by the netlib routine dgesv, and some- 
times also dgesvx) to invert the resulting system. A discussion of the struc- 
ture of the matrix serves to sharpen our analytic understanding of the linear 
systems we are dealing with. For 3d problems, however, speed and memory 
considerations make it impractical to form and solve the full matrix. For the 
3d model we consider in the conclusion and for envisioned future 3d work, we 
have used, and plan to continue with the Krylov-based method GMRES, in 
which only the specification of a matrix-vector multiply need be implemented. 
For a Krylov method preconditioning strategies are often necessary to avoid 
stagnation of the iterative solver. 

Despite the name of the method, and despite its intuitive appeal, it is not 
guaranteed that IPC, especially for rectangular domains, actually improves 
conditioning (with respect to the problem of matrix inversion). If this is an 
issue in two dimensions, it is sure to be even more problematic in three di- 
mensions. Since we cannot guarantee that our implementation of IPC will ac- 
tually improve conditioning, it would be more appropriate to call the method 
a "sparse formulation." To refer to the very specific IPC technique it will be 
convenient, however, for us to continue to use the term "preconditioning." 
For simplicity, and to have a uniform treatment, we will use the IPC method 
for all our subdomains in the 2d model. It should be understood that this 
uniformity is not necessary. We could treat the elliptic region via a different 
method (spectral or pseudospectral). In any case, as we verify with numerical 
experiments, the method yields impressive global accuracy in solving the 2d 
HRWE on our two center domain. 



3.1 Basic formulas for Chebyshev polynomials 



Although Ref. 40| considered general orthogonal polynomials, we work solely 
with Chebyshev polynomials, a classical system of orthogonal polynomials 
with many useful properties and applications [i^, H^. Here we collect only 
those properties relevant for our discussion, mostly following 40, 3, 45|. The 
degree-n Chebyshev polynomial Tn{^) is defined by Tn{^) = cos(?t, arccos(.^)) 
for —1 < ,^ < 1, showing that we may consistently set T_n{^) = T„(0- In 
our application ^ depends on one of the coordinates, say C,{r) = (2r — r^ax — 
'^min)/('"max — ^min) in an outer Spherical shell or annulus. The T„(0 ^.re solu- 
tions to a singular Sturm-Liouville problem, and therefore particularly suited 
for approximating solutions to differential equations on [—1, 1] with a wide 
class of boundary conditions. The T„(0 ^ire orthogonal on the interval [—1, 1] 
with respect to the weight function (1 — ^^)~^/^. 
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We have To(^) = 1, T'i(^) = ^, and the following identity: 

2eT„(0 = T„+i(0+T„_i(0. 

If we denote by T(0 the row vector [To(0, ^2(0, • 
be expressed as 

eT(0 = T(OA, 

wher^ 

^0 i ■■ 
1 i ■■ 
i I ■■ 
A= 00i0|0-- 
1 i ■■ 
i ■■ 



(23) 

then (23) may 
(24) 



(25) 



Note that our convention here is to label matrix rows and columns starting 
from rather than 1. Since A is tridiagonal, the matrix A*", representing left 
multiplication by ^""^ through the formula ^"^T(,^) = T(,^)yl"\ has bandwidth 
2m+l. Moreover, banded matrices p{A) similarly correspond to multiplication 
by any polynomial p{C)- The most important special case, is the matrix At^ = 
Tm{A) of bandwidth 2m + 1 corresponding to multiplication by the Chebyshev 
polynomial Tm{0- Fo'^ m = 1, the matrix A = At^ is given by Eg. (25). More 



generally the entries of At^ may be gathered from the identity |40l. |4 



2T^(0T„(0 = T„+„(0 +T|„,_„|(0 . 



(26) 



Next we consider the formal expansion 

oo 

n(0 = E^n^n(0, (27) 

n=0 

where, owing to the orthogonality of the Chebyshev basis, the expansion co- 
efficients have the analytic values 42 

Un = ^nj\0T40^y=^- (28) 

The prefactor /t„ is 1 /tt for = and 2/7r for n > . We denote by u the col- 
umn vector of expansion coefficients [uq, ui, U2, - ■ ■ f. In the vector space of ex- 



^ For now we proceed with infinite-size matrices. However, numerical approxima- 
tion necessarily entails suitable truncations. Our practice will be to use the same 
symbols for such truncations, and alert the reader whenever the viewpoint shifts. 
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pansion coefficients multiplication by ^ is effected through multiplication by A, 
that is to say = Au represents the expansion coefficients [mq, ^i; ^2; ^35 ' ' ' ]* 
for u^{^) = Similarly, in the space of expansion coefficients, multipli- 

cation by any polynomial p(^) may be represented via left multiplication by 
p{A). 



Another well-known formula [i^, 42, 44 



2(n+l) 2(n-l)' 



(29) 



yields another matrix identity 



T(0 = T'(Oi?[i] 



(30) 



in terms of an integration matrix 





1 -i 
Oi -i 
i -i 
i 



Jo 



(31) 



Here the subscript [1] merely emphasizes that the ffist row of consists only 
of zeros. Adopting a ffist row of zeros is a choice, permissible since Tq(^) = 0. 
Free, or open, rows of zeros within integration matrices is a central feature 
of the spectral method presented in 40|], and following that reference we will 
eventually exploit such freedom to enforce boundary conditions. In matrix 
form the formula 



+ 



4(n + l)(n + 2) 2(n2 - 1) A{n-l){n-2) 



(32) 



reads 



T(o = r'iOBl 



2]' 



(33) 
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where now 



Bl 







1 -i 

4 " 6 












J_ 

24 









I, 



h 



0^ ^ -h 



-h 



(34) 



Single and double integration of u then correspond to i?[i]U and -B^jU. 



Consider the matrices D and corresponding to differentiation in the Cheby- 
shev basis. The entries of these matrices stem from recursive use of (29) for 
D and (32) for D^. Although we will not require the precise forms of these 
matrices, we nevertheless list them here for completeness: 



D 



10 3 


5 





7 





4 8 





12 





16 


6 


10 





14 





8 





12 





16 





10 





14 











12 





16 











14 

















16 



(35) 
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4 32 108 256 

24 120 336 

48 192 480 

80 280 

120 384 

00000 168 

00000 224 

00000 

00000 



(36) 



The matrices D and obey 



(37) 



Notice that the first column of D and the first two columns of are (neces- 
sarily) zero. We might accordingly represent these matrices as D(^i) and to 
highlight these facts, but will not employ this notation. Whereas the matrices 
and B'^^ are banded and sparse, D and D"^ are clearly upper triangular 
and dense. To render the linear systems we encounter in banded form, we will 
exploit the identities, 4^ 



Bii]D = /[I], 



B[2]D = B[2]. 



(3J 



Here the notation J[i] means the identity matrix with the first row set to zero; 
I[2] means the first two rows are zero. We stress the crucial point that, for 
the ordering of matrix products in Eq. (38), the identities (38) also hold for 
(A^ + l)-by-(iV + 1) truncations of Bf^^, D, and D^. 

In order to enforce boundary conditions we introduce both Dirichlet and Neu- 
mann vectors (radiation boundary conditions require a bit more thought and 
will be put off until later). The Dirichlet vectors are 



5+ 

6" 



1,1,1,1... 



-1,1,-1,... 



To(i),ri(i),r2(i),T3(i),-- 

ro(-i),Ti(-i),r2(-i),r3(-i) 



(39) 
(40) 



while the Neumann vectors are 



V 



0,1,4,9, 



-1,4, 



T^(1),T;(1),T^(1),T^(1),-- 

r^(-i),T((-i),r^(-i),T^(-i), 



(41) 
(42) 



18 



The boundary conditions we encounter are easily expressed in terms of these 
vectors. For example, 5~-u = is a homogeneous Dirichlet boundary condition 
at the left endpoint, and i/"*" ■ u = is a homogeneous Neumann boundary 
condition at the right endpoint. 

We conclude this preliminary subsection by addressing a few practical mat- 
ters. First, our discussion so far assumes the standard interval [—1,1]. Mod- 
ifications, albeit trivial ones, of the matrices A, B[i], B^2] vectors z/^ are 
needed for different intervals. For example, if we use the T„(^(r)) to approx- 
imate functions on r G [rmin, '"max], then to represent double integration that 
corresponds to the physical interval we must send 0.25(rjnax— Tmin)^-B|]- 

Similar scalings must be made for A, B[i], and i^^. However, we shall retain the 
same symbols for matrices which are so scaled. Second, although we may the- 
oretically consider an infinite expansion (27) for a given function u (perhaps 
a solution to a differential equation), numerical applications entail a suitable 
truncation Vnu{C) = Yln=o''^riTn{0^ where Vn represents the truncation op- 
erator. Moreover, most applications do not work directly with the analytic 
expansion coefficients (28). Rather, one typically approximates the integral in 
(28) via the quadrature rule stemming from Chebyshev-Gauss-Lobatto nodes 



(see, for example, (4^). This process introduces aliasing errors, and results in 
discrete expansion coefficients u'^^'^'^'^^'^ . A more theoretical treatment would, 
when appropriate, draw a careful distinction between the analytic and discrete 
expansion coefficients, but this issue is not of primary importance to us. 

3.2 Sparse formulation on rectangular domains 

In terms of co-moving Cartesian coordinates, 

X = r cos (f = r cos{(f) — Qt) , y = r sin (f = r sin{(j) — Qt) , (43) 
let the scalar field ip{x,y) obey the inhomogeneous 2d HRWE 



dl + d^y-n^[xdy-yd. 



^ = g (44) 



on a rectangle R = [xmm, a^max] x [ymin, 2/max]- We have also worked with the 
shifted equation resulting from the mapping x^x — c, y—^y — d, where 
(c, d) is the center of the rectangle. For simplicity, we do not consider a shift 
here. Via repeated use of the Leibniz rule, we write (44) as 



diil - nY) + d'yil - n'x') - n\d,x + dyy - 2d,dyxy}\ ij = g, (45) 
in preparation for our sparse-matrix spectral approximation. 
At the theoretical level, we may represent our solution in terms of a double 



19 



Chebyshev expansion, 

^{^^y) = J2T. ^nMa^mUviy)), (46) 



oo oo 



n=0 m=0 



where iC{x)yV{y)) is a mapping of our rectangle R onto the unit rectangle 
[— 1, 1] X [—1, 1]. To obtain the system of equations we solve numerically, we 
consider the truncated series 

N M 

VNM'^i^^y) = J2J2 '^nmTni^{x))T^{ri{y)). (47) 

n=0 m=0 

We represent the finite collection of expansion coefficients as a 1-vector 
V' = (^00, ^01, ■ ■ ■ , i^oM, ^10, V'li, ■ ■ ■ , ^iM, ■ ■ ■ , i^m, ^Ni, ■ ■ ■ i^Nu) , (48) 



so that the components il^ik) = {'4')k are determined by the direct product 
representation 

V^(n(M+ 1) + m) = V^,^. (49) 

We then consider the approximation of (45) in terms of i/' and suitable trun- 
cations of spectral differentiation matrices. 



Dl ® (J, - fiMJ) + (4 - ^'AD ® Dl 



Vt^D^A^ ®Iy + h® DyAy - 2D^A^ (g) DyAy 



V' = g- (50) 



Here represents the (A^+l)-by-(A^+l) differentiation matrix in the Cheby- 
shev basis T„(^(x)), and A^ represents the (A^ + l)-by-(A^ + 1) matrix corre- 
sponding to multiplication by x, with similar statements for Dy and Ay which 
are both (M-(- l)-by-(M + 1). As mentioned at the end of the last subsection, 
to obtain these matrices appropriate scaling factors must be included with the 
straightforward truncations of the infinite-size matrices listed above. By the 
definition of the Kronecker direct product, we have, for example, that 

{Dl ® Al)in{M + 1) + m, kiM + 1) + p) = {Dl){n, k)iAl){m,p), (51) 

in the Matlab notation C{i,k) for entries of a matrix Cik- 

To achieve a sparse and banded representation of the approximation, we mul- 
tiply (50) by i3 = i?^[2] ®-^^[2] exploit the identities (38), thereby reaching 



+ -Sx[2] ® Byl2]Ay - 2B^12]A.J: ® Byl2]Ay 



V' = 5.[2]®5,[2]g- 

(52) 
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The coefficient matrix of the system is then 







^ 

" 24 





































^[2] 

6 





^[2] 
24 











^[2] 

16 





^[2] 

48 





48 





30 











I[2] 

80 





48 






+ 







AN{N-1) 



2(Ar2-i) 












■■ 


■ 











■■ 


■ 











■■ 


■ 





B[2] 





■■ 


■ 










■■ 


■ 













■ 



■ ■ ■ 5p2] 

where we have only exphcitly displayed the fi-independent contribution to 
the matrix stemming from the preconditioned Laplacian. This Laplacian con- 
tribution exhibits the main features of the overall matrix. Namely, that it is 
banded, sparse, and contains free rows of zeros. Each entry in either of the 
matrices above is itself an (M + l)-by-(Af + 1) matrix. Indeed, the J[2] and 
Bp] matrices are Iy[2] and i?^^] in the notation of Eq. (52), and in the y di- 
mension the truncation is m = 0, ■ ■ ■ , M. Overall, the coefficient matrix has 
(A^-|-l)^ such blocks. The ffist two block-rows in the matrices displayed above 
are empty, giving 2(M -|- 1) free rows of zeros. In each of the remaining — 1 
block-rows the ffist two rows are empty, giving us another 2(iV — 1) zero rows. 
Therefore, we have a total of 2{N + M) such zero rows at our disposal. 



In the free zero rows we insert the "r-conditions" 40|, that is, the boundary 
conditions. To illustrate, we here assume Dirichlet conditions '?/'(xmin, 2/) = 
^(xmax,2/) = /"^(l/), ^(a;,2/min) = h' {x) , ilj{x,y^^) = h+{x). These 
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boundary conditions can be approximated as 



M 



N 




(53) 



m=0 



n=0 



where one-dimensional Chebyshev projections appear on the right hand sides. 
There are 2(A^ + 1) + 2(M + 1) boundary conditions in (53), but we now show 
that they are not all linearly independent due to double counting of corner 
conditions. The value of ^/^(xmax, ?/maxl = ^"^(a^max)^ /"^(l/max) can be written 
as a linear combination either of the or of the . This implies a homoge- 
neous linear relationship between the summations in the first and the second 
set of equations in (53). There are three other such linear relationships that 
follow from the ways in which ^(a^max, l/min), V'(a^mm, ymax),^and 1p{Xmm,ymin) 
can each be expressed either in terms of sums of or of Therefore, the 
number of independent equations in (53) is 2(iV+l)+2(M+l)-4 = 2{N+M), 
precisely equal to the number of zero rows. 

All methods for eliminating the four linear dependencies in (53) should yield 
comparable accuracy for a numerical solution. However, dropping the four 
highest-mode equations {n = N for ± in the left equation, and m = M for ± 
in the right equation) is inconsistent, as confirmed by numerical experiments. 
Indeed, by throwing out the highest mode on all four edges, one loses infor- 
mation about the corner values. We have chosen to reduce (53) in a consistent 
way which is at the same time convenient for our direct product representa- 
tion of the rectangular region. We use all 2N + 2 of the first set of equations 
in (53), i.e. those equations involving the h^. The Dirichlet vectors 5^ associ- 
ated with these conditions are placed within the first two rows of each block 
in the coefficient matrix (which does not increase the bandwidth of the coef- 
ficient matrix beyond the block-diagonal). We must then dispense with four 
equations from the second set, two for — (left) and two for + (right), and so 
drop those equations for which m = M — 1 and m = M. This corresponds to 
dropping the two highest modes for both left and right Dirichlet conditions, 
although we have experimented with dropping other modes and, as expected, 
found little difference. From knowledge of all the (information determining 
all four physical corner values) and the for m = 0, ■ ■ ■ , M — 2 we may re- 
cover the four values for f^-i and /^j, using the four equations for the corner 
values expressed as linear combinations of the The Dirichlet vectors 6^ 
associated with the remaining left/right boundary conditions are then placed 
in the remaining 2M — 2 rows of the first two blocks (they reach across block 
columns and so affect the bandwidth of the coefficient matrix quite a bit). 
Values appearing on the righthand side of the above boundary conditions, all 
and /^^ except for and /f^, replace the appropriate zero entries of 

the source ^^^j O fij^jg. 

Although, our preconditioning of the rectangular operator has resulted in a 
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sparse matrix (while the original matrix corresponding to the Ds is dense), 
the issue of condition number is more subtle. Clearly, the Bs and Ds are 
in some sense inverses of each other. Notice that as an infinite dimensional 
matrix has the action: Qq~^ — > Qi , with denoting the vector subspace 
of spectral coefficients corresponding to Chebyshev expansion from degree k 
through degree p polynomials. Suppose we consider the (A^ + l)-by-(iV + 1) 
truncation of B^i] , and further that we delete the first row and last column of 
this square matrix, thereby obtaining an N-hj-N matrix B. Likewise, we take 
the (A^+l)-by-(iV+l) truncation of D, and delete its first column and last row 
to obtain D. Then we may view D : Q'^ Qo~^- Taken as square matrices 
with the same domain and range, B and D are nonsingular and inverses of 
each other, whence have the same condition number!^ Such an argument can 
produce nonsingular matrices B^ and D^, also inverses of each other with the 
same condition number. Therefore, although we are ignoring the critical issue 
of boundary conditions, passing from a coefficient matrix with symbolic form 
/ ® + ® I (corresponding to the Laplacian part of the operator) to 
one with the form B'^ ® I + I ® B"^ is not clearly advantageous insofar as 
the conditioning of the resulting linear system is concerned. Nevertheless, one 
might expect that a better distribution of eigenvalues for the form B"^ ® I + 
I ® B"^ would lead to faster convergence were we using an iterative solver such 



as GMRES 46 



3.3 Sparse formulation on annular domains 

In our 2d HRWE Lip = g, we now take x = a + pcosO and y = h + psinO, 
where the "hole" is located at (a, 6), which is either the center of annulus H 
or the center of A. In terms of 

F(^) = asin^-6cos^, G(^) = a cos ^ + fesin ^, (54) 

the HRWE operator is 

L = dl + p-'d, + p-^dl - [Fi9)d, + (1 + p-'Gi9))de] ' . (55) 

To prepare for the integration preconditioning, we use F'{6) = G{9) and 
G'{6) = —F{6) along with repeated appeals to the Leibniz rule in order to 
obtain 

p'Lij = p^g (56) 

with 



Recall that the condition number of a matrix A is k{A) = \\A ^||||A||. 
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(1 - fi'F^) + 9pp[-3 + fi'(2F^ + + pG)] 
+dl[l - n\p + + n^depF - 2n^dedppF{p + G) 
+1-Q^{G^ + pG). 



(57) 



We now represent the solution on the annular subdomain in terms of a trun- 
cated Fourier-Chebyshev expansion, 

N 

'PM,N'ip{p, 6) =J2 ^OnTni^ip)) 

+ E E [v^2fc-i,nCos(fce) + V2fc,nSin(A:^)]T„(e(p)), (58) 

k=l n=0 

where for simplicity we have here chosen M even. For the direct product 
representation, 

t/^((iV+l)m + n) =^^„, (59) 

which means each Fourier mode corresponds to its own (A^ + l)-by-(A^ + 1) 
Chebyshev block. We therefore have the following matrix representation C of 
p'L: 



£={Ie-n'f')®DlAl 

-Q^ (G - 2Def) ® DpAl - Q^D^g ® 



DoF- h + 2Di]G 



(S)Ap 



h + D'i) (le-n'G^)®Ip. (60) 



where san serif F and G denote matrices in the Fourier sin/cos basis corre- 
sponding to multiplication by F and G. The entries of these matrices are 
determined by standard trigonometric addition-of-angle formulas, such as 
2 sin a cos /5 = sin(Q; + j3) + sin(a — /3). For the scenario of the 3d HRWE, 
such a spectral approximation for a spherical shell around an inner hole would 
be rather more problematic (even if the analogous matrix is not explicitly 
formed). Indeed, for that scenario we would need to contend with Wigner- 
Clebsch-Gordon coefficients arising from products of spherical harmonics Y^m- 

Since the spectral differentiation matrices Dp and Dp are dense upper trian- 
gular, passage to a sparse-matrix formulation of the problem for an annulus 
necessarily requires that we apply the integration matrix -B^[2]- By contrast, 
since the matrices Dg and Dg are already banded or diagonal, we will achieve 
a sparse formulation whether or not we precondition with 6 integration. In 
fact we have employed 9 preconditioning, since it might well yield a better dis- 
tribution of eigenvalues and better equilibration properties. Nevertheless, for 
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simplicity will here ignore IPC in 6, which, in any case, is logically different 
from IPC in p in that no boundary condition is associated with the periodic 
direction. The preconditioning integration matrix is then B = Iq ® -B^j2] ' ^"^^ 
its application onto (60) yields 



EC 



Bp[2]Ap 



-Q^ (G - 2DeF) ® E^pjA^ 



2Dl) G 



^ Bl[2]Ap 



B 



P[2] • 



(61) 



The right hand side of 



p^g is now represented by B ■ {Iq ® A^)\ 



{Ig ® B'f^,2]A'l)g. For the outer annulus we have (a, 6) 



and 9 = (p. Also for this case F(9) = 



(0, 0), and so p = r 
G{9), and C in (60) reduces to 



L. 



(62) 



We therefore find 



BC = I^® Ir[2\Al - 3/^ ® Br[2]Ar + ® B^^^ (/, - fi^A^j +1^® E^^j (63) 

as the preconditioned matrix for the outer annulus. 

Turning to the issue of boundary conditions, we first note that for EC the first 
two rows of each block have all zero entries. Into these rows we therefore place 
the r-conditions, 

N 



n=0 



± 



mn^n 



k 



± 



(64) 



(and similarly for Neumann conditions). Here is the Fourier transform of 
the boundary conditions h^{9), for example with h~{9) = ip{pmin,9). Notice 
that the correct number of zero entries in the preconditioned source correspond 
to the inhomogeneity h^. For the radiation "p, q boundary conditions" (19) 
we are only concerned with the outer annulus and the simpler operator (63). 
In this case, the overall (iV + 1)(M + l)-by-(A^ + 1)(M + 1) matrix BC is 
block diagonal, with (M + 1) blocks. Each block has the structure 








BC 



(65) 



for some Fourier wave number k. Here the represents a row of zeros, and 
BC^ is a nonzero (A^ — l)-by-(A^ + 1) submatrix. As the radiation boundary 
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conditions couple sine and cosine modes of the same wave number, we must 
consider the two consecutive blocks (one cosine and the other sine) associated 
with Fourier wave number k ^ 0. As depicted in the last equation, each of 
the two blocks has all zero entries in its first two rows. The p, q boundary 
conditions are enforced by choosing the overall block neighborhood as follows, 



s~ 







Rv+ + q5+ 


BC^ 








5- 


Rv+ + q5+ 


-p5+ 





BC^ 



(66) 



where represents either a row or a (A^ — l)-by-(A^ + 1) submatrix of zeros. 
Boundary conditions for A; = (the zero mode) are easier to enforce (only a 
single block need be considered) and are handled similarly. 

3.4 Gluing of subdomains 

So far we have described individual rectangular and annular subdomains (and 
their associated r-conditions) as if these subdomains were decoupled. In fact, 
we "glue together" all or most of the subdomains shown in Fig. 1. This glu- 
ing takes two forms: (i) imposing matching conditions for adjacent rectangles 
and (ii) imposing matching conditions for the overlap between an annulus and 
a set of rectangles. Before describing each case in more detail, we comment 
on how such gluing is reflected in the overall linear system. Let and 
represent the vectors of Fourier-Chebyshev expansion coefficients associated 
with the spectral representation of the solution on the annuli H and A. Sim- 
ilarly, let if)^ represent the vector of double-Chebyshev expansion coefficients 
associated with the spectral representation of the solution on the jth rectan- 
gle (with 1 < j < 8). The overall set of unknowns is then the concatenation 
^ = (■0-'^, ■0-^, -0^, ■0^, i/j^, ■0^, ■0^, '0'^, '0^, ■^^)*, which satisfies (for the linear 
problem) the spectral matrix form of Eq. (1) 

= BQ, (67) 

where ^ is a similar concatenation of the sources g on the individual sub- 
domains and the B indicates integration preconditioning on all subdomains. 
Symbolically, the coefficient matrix Ai is BC, here with C standing for the 
spectral representation of the HRWE operator L on the whole two center 
domain. 
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Each of the ten subdomains (annuh A and H, as well as rectangles 1-8) in 
Fig. 1 are represented by one of ten super-blocks {H-H, 1-1, 2-2, ■ ■ ■ , 8-8, 
A-A) which sit along the diagonal of the overall super matrix A4 representing 
the PDE on the whole two center domain. We use the term "super-block" 
here since the matrix corresponding to each subdomain arises, as we have 
seen, from a direct product structure (and so could be viewed as already in a 
block form). The supplementary equations needed for gluing are placed within 
existing zero rows in the same manner as with the r-conditions. However, 
the gluing conditions stretch beyond the super-block diagonal, since they are 
linear relationships between the spectral expansion coefficients on two (or 
more) separate subdomains. For example, the gluing together of subdomains 
1 and 2 (which share a common edge) involves not only filling rows within 
the 1-1 and 2-2 super-blocks along the diagonal of Ai, but also filling rows 
within the 1-2 and 2-1 off-diagonal super-blocks. 



3.4-1 Gluing of rectangles to rectangles 

For rectangles which meet at an edge we require both continuity in and 
its first derivative dip/dv (normal to the matching edge). We impose these 
requirements strongly, that is to say at the level of the numerical solution itself. 
Consider, for example, rectangles 1 and 2 in Fig. 1, which as indicated share the 
edge y = Pmin, where pmin is the radius of the inner hole (the depicted excised 
region). We require that values of 4'{x,pmin) and dip / dy{x, p^ain) agree at the 
Chebyshev-Gauss-Lobatto collocation points x{S,i) whether these values are 
computed using the spectral coefficients of rectangle 1 or those of rectangle 2. 
For simplicity, we assume that the number A'^i of spectral elements for rectangle 
1 is the same as that N2 for rectangle 2, so that the matching equations are 
simply 

Ml ^ M2 _ M\ ^ M2 ^ 

^Im^m = i^lm^m, Y ^lim^^^^^ = " II ^lm(^2J^m^ (68) 

m=0 771=0 m=0 m=0 

for each value of n. The a factors here are the scalings of the Neumann vectors 
that are necessary since the range of y may not be [—1,1]. (See the discus- 
sion following Eq. (42).) These matching conditions are refiected in the overall 
matrix M. as follows. As the super-block corresponding to each of the subdo- 
mains 1 and 2 has been preconditioned in the described fashion, each has a 
collection of zero rows in which we place the matching conditions. In, say, the 
zero rows belonging to the super-block 1-1, we insert the first set of conditions 
given in (68). In the zero rows belonging to the super-block 2-2, we similarly 
place the Neumann conditions, the second set of conditions given in (68). We 
note that this filling of zero rows to achieve the required matching does not 
affect the inhomogeneity, as these are homogeneous conditions (a linear sum 
of expansion coefficients for one subdomain plus a linear sum of expansion 
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coefficients for another is set equal to zero). 

In practice, the gluing procedure is plagued by the same sort of issues that 
arose in connection with Eq. (53). Ultimately due to redundant counting at 
corners, the full set of matching conditions responsible for the gluing of all 
subdomains (including annuli) are not linearly independent. Thus we find 
that the number of apparent matching (and boundary) conditions is greater 
than the number of zero rows available in the super matrix Ai to hold the such 
conditions. As we did for a single rectangle, we have chosen to consistently 
eliminate linear dependencies at left-right edges/interfaces, as illustrated in 
the next paragraph. The number of zero rows in A4 always equal to the number 
of linearly independent matching (and boundary) conditions, and in principle 
any two ways of implementing them should be equivalent. While the choice of 
implementation does of course affect the structure of the super matrix Ai, it 
should not greatly affect the accuracy of our numerical solutions. 

Let us turn to our illustrative example, the matching, along the vertical edge 
X = xh + Pmin, of rectaugles 2 and 3 in Fig. 1. In the matrix for rectangle 2 
(super-block 2-2 of the overall super matrix Ai) we have reserved 2M2 — 2 
rows in the first two blocks for enforcing boundary (or matching) conditions 
at the left and right edges. We will now have M2 — 1 of these available for 
matching at the left edge of rectangle 2. Likewise, in the matrix for rectangle 3 
(super-block 3-3 of the overall super matrix Ai) we will have M3 — 1 available 
rows for the matching along the right edge of rectangle 3. We take Mi = M2 = 
M, so that we have a total of 2M — 2 rows in the overall super matrix Ai 
available to explicitly enforce the matching of rectangles 2 and 3. The full set 
of equations for gluing along this edge are analogous to (68), but with A^s and 
Ms interchanged, 

£ = ^Im^n, ^lm(^2l^n = - Y i^lm(^3^n ■ (69) 

n=0 n=0 n=0 n=0 

We follow a protocol similar to the convenient one outlined for fixing Dirichlet 
boundary conditions for a single isolated rectangle. Whereas we have chosen 
to devote a zero row in Ai to each of the 2A^ + 2 matching conditions (68) for a 
top-bottom gluing, for a left-right gluing we devote only M — 1 zero rows for 
each set in (69), that is 2M — 2 rows in all. This means dropping the last two 
(highest mode) equations from each set. Since we work with the transform of 
the matching conditions, dropping the highest two modes does not mean we 
are not enforcing matching at a particular point along the interface. Moreover, 
although we have chosen not to reflect all of the explicit matching conditions 
(69) in At, extra matching of rectangles 2 and 3 is afforded by other matching 
conditions in the overall linear system. For example, along both the top edge 
of rectangle 2 and top edge of rectangle 3, we enforce a full set of boundary (or 
matching) conditions with Dirichlet data coming from annulus A. Since that 
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data will be smooth along the full length of the combined top edge for both 
rectangles 2 and 3, these top-edge boundary conditions effectively yield extra 
matching for the 2-3 vertical gluing, in particular enforcing both continuity 
and differentiablity at the physical corner point common to each rectangle. 

Another complication in gluing is the implementation of symmetry. Notice 
that the left edge of rectangle 4, for example, must be flipped and then glued 
to the left edge of rectangle 7. Similar flipped gluing must also be considered 
for subdomains 5 and 6. We enforce such flipped gluing of edges by demanding 
continuity in ip only (since region 7 is already matched to region 6), and filling 
the available zeros rows of the 4-4 super-block of A4 accordingly (as in this 
case the 6-6 and 7-7 super-blocks would no longer have appropriate free rows 
available). 

Let us provide a brief sketch of rectangle gluing without the "conforming 
assumption." Were we gluing, say, rectangle 1 to rectangle 2 (a top-bottom 
glue) with Ni ^ N2, we would deal with more complicated equations than 
(68). The second equation in (68), for example, could be written as 

Ml _ 

IplmO'l^m = hn, (70) 

m=0 

where is simply the explicit sum on the righthand side when A'^i = 
However, were A'^i N2, it would then arise as the x-Chebyshev transform of 
the vector 

N2,M2 _ 

K= i: ^i,T,{e{xn)) 

p,q=0 

Here the the Ni + 1 (scaled) Chebyshev-Gauss-Lobatto points for rect- 

angle 1, and {^'^{x),'r]'^{y)) are the linear functions mapping rectangle 2 to the 
square [—1, 1] x [—1, 1]. As there would be A^i + 1 of the equations (70), they 
would be placed in zero rows corresponding to the 1-1 super-block of Ai, and 
would stretch across the 1-1 and 1-2 super-blocks. Equations for continuity 
similar to the first set of equations in (68), but using the x-Chebyshev trans- 
form for rectangle 2, would determine entries in Ai stretching across the 2-1 
and 2-2 super-blocks. The roles of 1 and 2 could be interchanged. The situ- 
ation is nearly the same for a left-right gluing, modulo the issue of dropping 
the two highest modes. 

The rows in the overall super-matrix Ai responsible for all of the rectangle-to- 
rectangle gluing result in Ai having a large condition number K{Ai). Indeed, 
for the truncations we later consider in our numerical experiments, the matrix 
for the elliptic region (which includes the H H super-block as well) has a 
reciprocal condition number RCOND (a diagnostic determined by dgesvx) in the 
10~^ to 10~^° range. As we will see in Sec. 4.2, such a large condition number 
leads to an accuracy problem for higher truncations. In an attempt to mitigate 



(71) 



11=0.-. 
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this problem, we have experimented with other strategies for rectangle gluing, 
for example, imposing point-space rather spectral versions of the matching and 
boundary conditions. Another possibility, suggested by a referee, and for which 
we had earlier carried out some numerical experiments, is to allow for overlap 
between adjacent rectangles. In this case the gluing procedure is similar to the 
described rectangle-annulus gluing, with the key advantage that all Neumann 
vectors can be done away with (this could lead to improved conditioning). 
The disadvantage is that significant overlap is typically required for good 
conditioning, and, therefore, more of the physical space is doubly covered. 
Beyond an increase in resources, such double cover also further complicates the 
gluing of the overall concatenation of rectangles to the annuli. Nevertheless, 
this idea deserves further study. To sum up, we have not exhaustively tried 
every approach, and our current one is quite possibly not optimal, although 
we believe it will be hard to do significantly better. Nevertheless, by switching 
from dgesv to dgesvx, we are able to achieve the requisite accuracy for our 
experiments. The routine dgesvx includes equilibration (which changes RCOND 
to the 10~^ to 10~® range) and iterative refinement ji^t] . Moreover, with our 
approach we have carried out the experiment ( "L = 5 with linear mappings" ) 
documented in Sec. 4.1 of [l^, the Laplace equation for the exact solution 
ln(x^ + y^) on a 10 X 10 square with an excised 1x1 inner square "hole." With 
dgesv we attain a best error measure which is about three orders of magnitude 
larger than the corresponding one reported in Fig. 3 of that reference. However, 
with dgesvx we find nearly identical results, and even a slightly better best 
error measure. 



3.4-2 Gluing of annuli to rectangles 

The annuli H and A depicted in Fig. 1 overlap multiple rectangles, and for 
this overlap the issue of gluing is complicated. Since the issue is essentially the 
same for the gluing of H to rectangles 1-8 or A to rectangles 1-4 and 6-8, let 
us here focus on the first case. For example, part of the outer edge of annulus 
H sits in rectangle 1; we require that the values of ip along this outer edge of H 
agree whether they are computed with the spectral representation oiip on H 
or with the spectral representation of il^ on rectangle 1. Similarly, the values of 
ip at the left edge of rectangle 1 must be the same whether computed using the 
rectangle 1 or if spectral representation. That is to say, along these boundaries 
we demand agreement in the numerical solution whether determined by the 
expansion coefficients i/?^ of the annulus or the coefficients for rectangle 1. 
An essential difference should be noted between this kind of gluing, and the 
gluing of two rectangles at an edge: Here we use Dirichlet matching along two 
curves rather than Dirichlet plus Neumann matching along a single shared 
boundary curve. 

For the outer circular boundary p = pmax of the annulus H, we enforce match- 
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ing at each of a set of Fourier collocation points {9m '■ m = 1,2, ■■ • ,Mh} 
through the following set of equations: 



Y: V'LC = hi, (72) 

n=0 



where now the h'^ are not fixed Dirichlet values. They now arise via Fourier 
transform of 

hl= ^PWe{Xm))W{ym)). (73) 

p,q=0 

Here the jth rectangle contains the point {xm,ym) corresponding to the col- 
location point {pma.x,9m)- These equations are placed within the zero rows of 
the H H super-block. The explicit matching equations can be captured by 
expressing the Fourier transform /i+ = J2k=i ^mkh^ as a matrix relation in 
terms of the ipj^^. 

The matching along one of the inner edges of a rectangle follows a similar 
pattern. As a concrete example, take again the inner edge of rectangle 1, and 
the condition 

E V'LC = fm, (74) 
n=0 

where the fm are now not fixed Dirichlet values. Rather, they arise as the 
y-Chebyshev transform of the vector 

Nh _ 

fm = E^S^p(^(P™)) 
p=0 

\Mh Nh ^ 

+ E E K-l.pCOs(fc^„0 + ^l^^MkOmpMprn))- (75) 
k=l p=0 



Here the point {pm, Om) corresponds to a Chebyshev-Gauss-Lobatto collo- 
cation point (a^min, l/('7m)) along the inner edge, and we have chosen Mh 
as an even integer for simplicity. (Note that the values of 9m here, with 
< m < Ml have no direct relationship to the Fourier collocation points 9m 
used in Eqs. (72) and (73).) Again, via a matrix representation fm = J^eCmifi 
of the transform, we may express this matching condition more directly. In 
any case, these conditions are placed within the available zero rows of the 1-1 
super-block of A4. 
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3.5 Nonlinear model 



We turn to the 2d nonlinear HRWE discussed in ISlI, that is 



Lij + 7ih{^)=g, h{ij) = -^^^. (76) 



Certain aspects of our treatment of the nonhnear term will likely not gen- 
eralize to 3d. Indeed, we have retained a pure s pec tral method (coefficients 



as unknowns), whereas most modern approaches (46| towards solving nonlin- 
ear equations via Newton-Krylov spectral methods have centered around the 
construction of preconditioners within the context of pseudospectral methods 
(point values as unknowns). Nevertheless, as discussed in the concluding sec- 
tion, we expect that some aspects of our approach will carry over to 3d insofar 
as an outer spherical shell is concerned. 

Again, we let ^' represent the vector of unknowns, that is the overall con- 
catenation of the spectral coefficients on all subdomains. Using the inverse 
spectral transformation available on each subdomain, we may produce from 
a collection of point-space values ^ = JF~^^' defined on each subdomain's 
spectral grid. Here we are using JF^^ to formally denote the process of inverse 
transformation on all subdomains. Our problem now is to solve 

W[^= + r]Bn(j^~^^) - Bg = 0, (77) 



here with B representing the application of integration preconditioning on all 
subdomains, Q the concatenation of all source coefficients g, and Ai = BC 
the coefficient matrix for the linear HRWE on the glued-together two center 
domain. The term jF/ifjF"^^) = J^h('^) is the spectral representation of 



the nonlinearity, where ^(^) is constructed pointwise on each subdomain's 
spectral grid. 

Solving this nonlinear system of equations requires a Newton-Raphson itera- 
tion, and we briefly describe how this has been implemented. First, we note 
that by using both (26) and addition of angle formulas, such as 2 sin a cos P = 
sm{a + (3) +sin(a — /9), we may convert the spectral coefficients u representing 
a function u on a particular subdomain into a matrix C„ whose action on the 
subdomain's vector space of spectral coefficients corresponds to multiplication 
by u in physical space. For example, on a rectangular domain, the resulting 
matrix would arise as 

N M N M 

E E UnmTniax))Tm{viy)) ^ = E E ^nmAr^ ® , (78) 

n=0 m=0 71=0 m=0 

where At„ = T^A,) is {N + l)-by-(A^ + 1), At^ = Tm{Ay) is (M + l)-by- 
(M -|- 1), and each matrix has been appropriately scaled (as discussed at the 
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end of Section 3.1). For both annuh and rectangles we have written subroutines 
for fast computation of such matrices from input coefficients. Then, given a 
concatenation U = (u^, U"^, u^, u^, u^, u^, u^, u'', u^, u"^)* of the spectral coef- 
ficients on all subdomains, we may also form a (super-block diagonal) matrix 
Cu whose action on the vector space of spectral coefficients for the entire two 
center domain corresponds to multiplication by u in physical space over the 
entire two center domain. 

Let represent the current Newton iterate. Then computing the next iter- 
ate = — Si^/ requires that we solve the linear equation 

{M + r/i3C/,,(*oid)} 5* = l^[*°'<^], (79) 

where the derivative function is 

= ^m±^. (80) 



The matrix C/i/(^oid) is constructed as follows: (i) obtain the point values ^"''^ = 
j^-i^oid^ (ii) construct the set of point values h'{'^"^'^), (iii) construct the 
concatenation of spectral coefficients U = J^h'{'^°^'^), (iv) using the coefficients 
U, form the matrix Cft/(^oid) as outlined above. We have taken advantage of 
the sparse nature of the preconditioner B to achieve fast matrix-matrix and 
matrix-vector multiplies, with the former required for quick computation of 
the matrix-matrix product i3C/j/(^oid) that is needed for each Newton iterate. 
Each such linear solve has been performed directly with dgesv. 

We have found that the described implementation of Newton-Raphson itera- 
tion works well for the 2d model, and for a wide range of parameter choices 
we have not needed to include line searches. Nevertheless, we recognize that 
such an implementation is not likely a viable option for 3d. Indeed, for spher- 
ical shells the scalar spherical harmonics Y^m are required for the angular 
basis functions. Beyond the memory burden of forming matrices for the 3d 
problem, even representation of multiplication by a nonlinear function in the 
spectral Y^^ basis is not as practical, because it would require extensive use 
of Wigner-Clebsch-Gordan coefficients (the simple multiplicative properties 
embodied in addition-of-angle formulas are no longer at one's disposal). How- 
ever, in a Krylov approach one might use the spherical harmonic transform 
and inverse transform to the same effect, carrying out all multiplications in 
point space. 
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4 Numerical analysis and results 



This section considers tlie conditioning of our numerical problem on the outer 
annulus A. It also presents our results in solving the HRWE on various sub- 
domains as well as on the whole two center domain. In the first subsection 
we consider only a single Fourier mode in the outer annulus, and we study 
the condition number of the matrix that must be inverted. This analysis also 
pertains to the analogous outer shell problem associated with solving the 3d 
HRWE. Indeed, as we have shown in the second section, the ODE and bound- 
ary conditions governing the 2d and 3d mode equations are essentially the 
same. The second subsection describes methodology and collects error tables 
associated with our numerical experiments in solving the linear and nonlinear 
2d HRWE. A third subsection presents preliminary results in solving the 3d 
HRWE. 



4-1 Condition numbers 



Consider the radial operator which stems from the linear helical operator via 
a Fourier-series transform; consider also the radial operator that arises in this 
way from the Laplace-Poisson operator. With k taken as the wave number, 
these radial operators are respectively (see left hand side of Eq. (5)) 

d^ d 

L'' = r^— + r A;2(l-fiV), (Helical operator), (81) 

dr^ dr 

with Dirichlet boundary conditions at the inner radius and exact outgoing 
radiation boundary conditions at the outer radius (D-R conditions), and 

d^ d 

P'' = r^— + r- fc^ (Poisson operator, = 0). (82) 

dr^ dr 

For the Poisson operator we examine both Dirichlet-Dirichlet (D-D) and 
Dirichlet-Neumann (D-N) boundary conditions. As the Laplacian is trans- 
lation invariant, the radial Poisson-Laplace operator can be put into the 
same Fourier form wherever the annulus lies. The HRWE mixes Fourier modes 
if the center of the angular coordinates is offset from the center of rotation, 
as in the case of the hole subdomain H. 



Here, however, we will work with the outer annulus A for which the HRWE 
itself does not mix modes, although, as we have seen, the radiation boundary 
condition does mix cos and sin modes of the same wave number. Therefore, we 
essentially have an ODE setting in which the conditioning properties of IPC 
are documented 40]. We point out that Ref. 4^] also refers to the conditioning 
of a "helical operator." However, the helical operator of Ref. i^l arises in the 
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context of the Navier-Stokes equations and is of course different from the one 
we consider. 



Before numerically probing the issue of conditioning for the operators above, 
let us present a heuristic argument as to why IPC achieves improved condi- 
tioning in certain ODE settings, including as special cases those of immediate 
interest to us. Suppose we have the equation 



where the Si(r) and so(r) are polynomials and the equation must be sup- 
plemented with appropriate boundary conditions. Here, one should think of 
^{r) as a mode, either ipk{i^) or iptmii^)- Performing the IPC technique on the 
ODE (83), we pass from a spectral approximation of the equation based on 
differentiation matrices. 



with two available rows of zero in which to put r-conditions specifying the 
boundary conditions to be applied. In either representation, Si{A) and 5*0(^4) 
are the matrices which correspond to multiplication by si(r) and so{r) in the 
Chebyshev basis. In the passage from (84) to (85), the structure of the co- 
efficient matrix is markedly improved. Indeed, in (84) the coefficient matrix 
features poorly conditioned dense matrices, save for So{A) which is sparse and 
banded. Moreover, the derivative operators are unbounded as grows. How- 
ever, in (85) the coefficient matrix is a perturbation of the identity featuring 
banded matrices, and the integration matrices remain bounded as grows. 

Let represent the matrix — built from Chebyshev differentiation matri- 
ces — which represents L'', ECJ' the corresponding preconditioned matrix, and 
likewise for and BV^ . Boundary conditions must be inserted into these 
matrices appropriately, and we will return to this point shortly. Our goal 
is to compute and compare condition numbers for the various operators and 
boundary conditions above. Given a matrix A, we shall work with the standard 
condition number k{A) = || A||2|| A~^||2 belonging to the matrix 2-norm. 

Recall that for the radiation p, q boundary conditions we must consider con- 
secutive blocks, since sin and cos modes of the same wave number are coupled. 
Suppose we consider the submatrix associated with the two consecutive blocks 
associated with Fourier wave number k. The unpreconditioned matrix for the 





D^ + DS^iA) + So{A) ^ = g. 




to the following matrix representation: 
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TruncSjtion 


Preconditioned 


Unpreconditioned 


o 


1 9 ^S'?Q n 9 fi^ 

iz.oooy ^^iz.Dj 


Q vT/ino-Ln"? 77/io_i_n'^'\ 
o./ ( "iue-i-uo \^o./ ( "ie-i-uoj 


iO 


9n fii^fiQ ('9n 
zu.Dooy (^zu.Dj 


o.ooz / e-t-u4i: (^o.oooe-f-u4i:j 








64 


51.1595 (51.2) 


1.3140e+07 (1.314e+07) 


128 


86.9027 (90.0) 


2.0748e+08 (2.075e+08) 


256 


156.0339 (156.) 


3.2976e+09 (3.298e+09) 



Table 2 

2-NORM CONDITION NUMBERS. D-D Poisson operator on [1,3] with k = 3. The 
numbers in parenthesis are those Hsted in Tables 1 and 4 of 40(] for the same 
problem. That reference's ho condition number is the 2-norm condition number. 

helical operator with D-R boundary conditions is then 













p5+ 


Ru+ + q5+ 












Ru+ + q6+ 


-p6+ 



with the understanding that the boundary conditions have been written over 
the last two rows of C'' in each block. However, as we have seen in (66), the 
preconditioned operator is 










Rv+ + q5+ 










5~ 




-p5+ 





EC'' 



now with the understanding that the boundary conditions have been written 
over the first two rows of each block. 



As an example of the various other matrices stemming from the Poisson oper- 
ator, consider the preconditioned matrix with Dirichlet-Neumann boundary 
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TrunccLtioii 


r^oisson [^u-L) ) 


i olsson \^u-iN j 




OZ 


1 '^7 9/1 f^9 


Q 07/1 p^o_un*^ 

o.U / 406-t-Uo 


y.oo4i:0e-t-u4 








9 7fizL9f3-l-nR 

Z. / U^Zcn^UU 


128 


159.5918 


6.3311e+06 


8.4756e+07 


256 


164.6046 


1.9820e+08 


2.6537e+09 


512 


180.0820 


6.2727e+09 


8.3984e+10 



Truncation 


Poisson (D-D) 


Poisson (D-N) 


Helical (D-R) 


32 


8.8450e+05 


6.3270e+06 


8.6394e+05 


64 


1.3791e+07 


9.8651e+07 


1.3471e+07 


128 


2.1776e+08 


1.5577e+09 


2.1270e+08 


256 


3.4610e+09 


2.4757e+10 


3.3805e+09 


512 


5.5191e+10 


3.9479e+ll 


5.3907e+10 



Table 3 

2-NORM CONDITION NUMBERS. Preconditioned (top) and unpreconditioned (bot- 
tom) operators on [5,15] for k = 2 and = 0.1. 

conditions, 



6- 





u+ 





BV'' 








6~ 





u+ 





BV^ 



Here we consider this decoupled direct sum of blocks only to make direct 
comparison with the matrices above. The condition number of the last matrix 
is of course the same as the condition number of either of its blocks. 

For our first investigation of condition numbers we take the Poisson operator 
associated with D-D boundary conditions. We assume /c = 3, a domain [1, 3], 
and consider truncations = 8, 16, ■ ■ ■ , 256, where each decoupled block is 
(A^-|-l)-by-(A^-|-l). We have computed the conditions numbers for these trun- 
cations in Matlab, and recorded our results in Table 2. Owing to the bound- 
ary conditions which are inserted into the zeroed rows of the preconditioned 
operator, we expect growth in the 2-norm condition number for our spectral 
matrix representation of the preconditioned operator as the truncation size 
N increases. The issue is significant for boundary conditions involving a Neu- 



37 



TrunccLtioii 


r^oisson [^u-L) ) 


i olsson \^u-iN j 




OZ 




i.uyzze-hUD 


i.yyuoe-t-uo 








Z.UUOOcn^UfJ 


128 


1.4650e+04 


1.0949e+06 


3.2270e+05 


256 


1.4673e+04 


1.1884e+06 


3.0785e+06 


512 


1.4731e+04 


4.2550e+06 


9.0043e+07 



Truncation 


Poisson (D-D) 


Poisson (D-N) 


Helical (D-R) 


32 


3.9154e+05 


2.9145e+07 


3.4677e+05 


64 


6.1040e+06 


4.5437e+08 


5.4046e+06 


128 


9.6378e+07 


7.1742e+09 


8.5335e+07 


256 


1.5318e+09 


1.1402e+ll 


1.3562e+09 


512 


2.4426e+10 


1.8182e+12 


2.1627e+10 



Table 4 

2-NORM CONDITION NUMBERS. Preconditioned (top) and unpreconditioned (bot- 
tom) operators on [5,150] for A; = 2 and $7 = 0.1. 

mann row vector, such as i/"*", whose entries grow hke with column location 
j. For k = 2 and Q = 0.1, Table 3 lists preconditioned (top) and unprecondi- 
tioned (bottom) condition numbers for the short domain [5, 15]. Table 4 lists 
the same numbers for the larger and more realistic domain [5,150]. For the 
helical operator the short-domain table shows that our preconditioning is only 
advantageous initially. As the truncation size increases, the r-conditions 
appear to have a dominant effect. (Since they mix the two blocks, for the 
helical operator such conditions would seem especially influential.) We loosely 
observe that the IPC technique only directly concerns the conditioning of that 
part of the coefficient matrix which stems from the operator itself, and not the 
boundary conditions per se. Over the whole range of truncations the large- 
domain table indicates that the condition number of the coefficient matrix is 
chiefly determined by the operator rather than the r-conditions. For this large 
domain we do discern the advantage of preconditioning the helical operator, 
although not as pronounced as for the Poisson operator. As a way of obtaining 
a more drastic improvement in conditioning of the helical operator, one might 
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explore either suitable equilibration tailored to the hr norms discussed in 
or enforcement of the radiation boundary conditions via a penalty method 44 



In Figures 3 and 4 we plot the singular values corresponding to the Helical (D- 
R) columns in Tables 3 and 4, respectively. In Figure 4, say, the top plot depicts 
the singular value distributions corresponding to the preconditioned matrices, 
for which condition numbers have been listed in the top leftmost column of 
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Fig. 3. Singular values for last columns of Table 3. For operators on [5, 15] 
with k = 2,n = 0.1. 

Table 4. The bottom plot shows the same data for unpreconditioned matrices, 
for which condition numbers have been listed in the bottom leftmost column 
of Table 4. In either plot, each curve corresponds to one entry of a column, 
where one should note that the number of plotted singular values is twice the 
listed truncation (each matrix consists of two modes with one block for each). 
Therefore, the longest curves corresponds to the 512 truncation, and so forth. 
Notice that the singular value distributions for preconditioned matrices are 
more clustered. 



4-2 Numerical results for the linear and nonlinear 2d HRWE 



Throughout this section we examine the 2d model, and, unless explicitly stated 
otherwise, all matrix inversions have been carried out with dgesv. Those in- 
versions carried out with dgesvx will be clearly indicated. 

Consider the linear equation 

= g ^(^-^^) [^(y, _ ^) _ ^(^)]_ (89) 
Xh 

where L is the HRWE operator and the inhomogeneity is comprised of two 
equal but opposite strength 5-functions, one located at (x, y) = {xh, 0) and the 
other at {x, y) = {—xh, 0). Assuming that exact outgoing radiation conditions 
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Fig. 4. Singular values for last columns of Table 4. For operators on 
[5,150] with k = 2,n = 0.l. 

are placed a.t r = R = npax? we may express the exact solution to (89) as the 
following infinite series [8|: 

ip{r,ip) = — y^^QJmjmQr^) Nm{rnQry)cos{m(p) + Jrn{rnQr^)sm{'mip) (90) 

m 

where the sum is over m = 1, 3, 5, ■ ■ ■ and r> is the larger of xh and r and 
similarly for r^. In what follows we always take xh = 2, Q = 0.1, and Q = 1. 
Although we shall not carefully study the pointwise convergence of the series 
(90), we note that it is poorly convergent whenever r is close to xh- In using 
(90) to generate approximations to field values ip for the experiments below, 
we have at times used 30 digit precision in Maple and included thousands of 
terms in order to sum the series with double precision accuracy. 

In order to make a direct comparison between the series (90) and numerical so- 
lutions, we now consider the following three linear and homogeneous problems. 



i. Elliptic problem on inner annulus. Let p represent the radial coordinate as 
measured from the source point (x//,0), so that the inner annulus corre- 
sponds to pmin < P < Pmax- We cousidcr the homogeneous problem Lip = 0, 
using the series (90) to place Dirichlet boundary conditions at p = pmin and 

P Pmax' 

ii. Mixed problem on outer annulus. We consider Lip = on the outer annulus 
corresponding to rmin = s < r < R = rmax, using the series (90) to place 
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9 

Fig. 5. CONSTANT-p SECTIONS OF ■ip{p,9). The figure depicts the solution (90) 
expressed in terms of the local polar coordinates (p, 6) about the hole. The top curve 
corresponds to pmin = 1> the middle to pmid = 1-5, and the bottom to Pmax = 2. 

a Dirichlet boundary at r = e and adopting exact outgoing conditions at 
r = R. In this setting, we assume e > xh (where e need not be small). 
Again, the light circle lies between e and R. 
iii. Mixed problem on two center domain. In this setting we consider the linear 
problem Lip = on the two center domain depicted in Fig. 1. We use the 
series (90) to place Dirichlet boundary conditions at p = Pmin (for both 
inner "holes" ) and adopt exact outgoing conditions at r = i?. 

In each case, we will generate a numerical solution, and then compare it with 
the series via various choices of error measure. Our numerical experiments for 
cases (i) and (ii) are meant to facilitate better parameter choices for scenario 
(iii) experiments. 

4-2.1 Linear elliptic problem on inner annulus. 

For the numerical experiment in (i) above, we use the series to generate Dirich- 
let data at pmin = 1 and pmax = 2. With a double Chebyshev-Fourier expansion 
and the described integration preconditioning, we numerically solve Lip = 



i 0.05 

n 
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20, 64 


3.192e-16 


1.399e-15 


1.352e-16 


7.766e-16 



Table 5 

Errors for the inner annulus problem. 

on this inner annulus. From the spectral solution we then generate numerical 
values for ip{pj^id,0), where pmia = 1.5 and 9 is sampled on a fine grid. This 
Pmid profile is then compared with the analogous profile given by the exact 
series (90). All profiles are depicted in Fig. 5, and Table 5 lists the associated 
errors for the pmid profile, where Nh denotes the number of radial Chebyshev 
elements, and Mh the number of angular Fourier elements. 

In this table, as in all remaining tables of the paper, the meaning of each of the 
error measures is as follows. The absolute sup error is simply the magnitude of 
the largest error found in the computed solution. For our spectral solutions this 
means that the solutions must be computed and compared on an evaluation 
grid. For Table 5 the grid was taken to be 1024 evenly spaced points of 6, at 
Pmid- The relative sup error is the absolute sup error divided by the maximum 
value found on the evaluation grid. The absolute rms error is found by taking 
the sum of the squares of all absolute errors on the evaluation grid, then 
dividing that sum by the number of terms that have been added, and finally by 
taking the square root. The relative rms error is the absolute rms error divided 
by the rms value of the solution on the evaluation grid. To evaluate the exact 
series on the evaluation grid, we have first evaluated it to high accuracy with 
Maple on a (uniformly spaced) Fourier grid of 128 points, subsequently using 
the corresponding Fourier interpolation onto the finer grid. Direct evaluation 
of the series on the 1024 point grid proved impractical. 

4-2.2 Linear mixed problem on outer annulus 

For this problem, in the region A of Fig. 1, we choose e = 4.5 and either R = 50 
OT R = 150. Since Q = 0.1, the light circle always lies within the corresponding 
annulus. Using (90) to generate Dirichlet data on the circle r = e is not 
problematic, since 4.5 > 2 and so the series converges reasonably well. We 
have used the FORTRAN routines RJBESL and RYBESL from net lib to build up 
the series (and also Matlab's Bessel routines for error analysis). In Table 6 
we list error measures in the outer annulus for various choices of Na (number 
of radial Chebyshev elements) and Ma (number of Fourier elements) . In each 
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32, 11 
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9.853e-16 
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Table 6 

Errors for the outer annulus problem. The top table corresponds to = 50 
and the bottom io R = 150. These numbers have been extracted from a 1024 x 
1024 uniform grid. 

case, for a fixed choice of Ma we have experimented to (roughly) find the 
smallest value Na which achieves the best possible errors. Figure 6 depicts the 
numerical solution %l}{r,ip) and the error A'?/;(r, (^), the difference between the 
numerical and series solutions, corresponding to the bottom table's second to 
last line. Notice that most of the error is concentrated near the inner boundary 
r = e, and this region determines the supremum errors in the tables. For a 
fixed particular number of angular Fourier elements, these values are, save for 
the last line, the same for both tables. 

We remark that owing to the quadratic growth inherent in the our radiation 
boundary conditions (based on Neumann vectors v^), 80 or 200 radial elements 
is likely excessive. We should prefer to further divide A into sub-annuli all 
glued together, with fewer Chebyshev elements taken on each sub-annulus. In 
future 3d work with numerical shells we plan to explore this possibility. 

4-2.3 Linear mixed problem on two center domain: iteration between the el- 
liptic region and outer annulus 

This numerical experiment tests the full multidomain spectral code. We use 
the series (90) summed in extended precision with Maple to place Dirich- 
let conditions on the inner holes (in practice, owing to the symmetry of the 
problem, on one inner hole H). We place radiation conditions at the outer 
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Fig. 6. Outer annulus problem. The plots depict the numerical solution and 
error for Na = 124 and Ma = 51. 

boundary r = R oi the outer annulus. We take the hole H and all rectangles 
1-8 of Fig. 1 as one glued domain, confined within a large square with side 
2L = 10. The outer annulus A is treated as a separate domain, and we iter- 
ate between solves on the hole-plus-rectangles domain and the outer annulus. 
More precisely, we first solve the hole-plus-rectangles configuration using the 
series (90) to place Dirichlet conditions on the inner boundary p = p^am of 
the hole and adopting Dirichlet zero conditions on the free edges of rectangles 
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Table 7 

Errors for the two center domain problem. These error measures have been 
extracted from the outer annulus only (again with a 1024 x 1024 uniformly spaced 
grid), although the numerical solution is generated on the entire two center domain. 
The last line of this table has been obtained with dgesvx (see text). Compare this 
table with the top table in Table 6. 

which overlap the outer annulus A. This multidomain numerical solution then 
provides inner Dirichlet conditions at r = e = 4.5 for the solve on the outer 
annulus. The solution for the outer annulus, in turn, provides a numerical 
solution from which we obtain corrected Dirichlet conditions for the relevant 
free edges of rectangles. While this iteration between solves on the hole-plus- 
rectangles domain and outer annulus is under way, we monitor the rms error 
between the inner profile '4){e, ip) from the outer annulus at the current and 
previous iterate. The iteration is stopped once this inner boundary rms error 
ceases to decrease, typically after 100 or so iterations. We could of course glue 
the outer annulus to the hole-plus-rectangles domain, and indeed we do so 
later when examining the nonlinear model problem, but here we use the de- 
scribed iteration, since we wish to demonstrate its stability. Table 7 lists error 
measures associated with this experiment with the outer boundary taken as 
R = 50. In the table the choices for Nh, Mh and Na, Ma stem from our ear- 
lier experiments on annulus domains. For simplicity, we have chosen the same 
numbers A^, M for the elements associated with the double Chebyshev expan- 
sion on each of the rectangular domains. For the table's last line, we have used 
dgesvx in lieu of dgesv, for the hole-plus-rectangles domain only. Without 
the extra accuracy afforded by dgesvx for the elliptic region, the numbers in 
the last line would be only marginally better than those in the next to last. 
See the final paragraph in Sec. 3.4.1 for further germane comments. 

4-2.4 Nonlinear mixed problem on two center domain 

For our final experiment in 2d, we examine the homogeneous [g = 0) nonlinear 
equation (76). For nonlinear solves we have always chosen to work with the 
fully glued two center domain {H + all rectangles 1-8 -t- A). Our numerical 
experiment with the nonlinear model is to repeat the mixed boundary value 
problem carried out for the linear model on the full two center domain. In 
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Fig. 7. Comparison of solutions. The plots depicts numerical solutions for the 
linear (top) and nonlinear (bottom) models. In each case the solution is not known 
in the immediate neighborhood of the holes, whence on these neighborhoods we 
have set the solution to a constant value in order to achieve nice plots. 

particular, we use the same Dirichlet conditions, from (90), at p = Pmin- A 
major difference is that now we solve simultaneously for the unknowns in all 
sub domains rather than by back and forth iteration with the outer annulus 
as we did for the linear problem. Since there is no exact solution available for 
comparison, we simply compare solutions with different resolutions. As before, 
the comparison is carried out using output associated with the outer annulus 
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Table 8 

Errors for the nonlinear two center domain problem. As with the hnear 
model experiment, these error measures are extracted from the outer annulus only. 
Error measures for the first two lines are taken with respect to the numerical solution 
corresponding to the next line. Error measures in the last line are taken with respect 
to the numerical solution which corresponds to incrementing all parameters by 2. 

on a 1024 x 1024 uniform grid. We take R = 50, rj = —10, and tpo = 0.25. 
With a numerical solution to the nonlinear problem, we may naively compute 
its "dipole strength" via Fourier transform of the solution restricted to the 
outer boundary r = R. The nonlinear term corresponding to our parameter 
choices results in a dipole strength which differs from dipole strength for the 
linear model {t] = 0) by about 18 percent, and this difference measures the 
nonlinear term's strength. Figure 7 depicts solutions for both the linear and 
nonlinear models in the vicinity of the holes, showing that even to the eye 
these solutions differ. Error measures are collected in Table 8. 



4-3 Preliminary results for the linear 3d HRWE 

We briefly describe two numerical experiments involving the 3d HRWE and 
meant to suggest that a multidomain spectral approach is viable for 3d PSW 
approximations, at least insofar as the post-Minkowski scenario is concerned 
and possibly beyond. The first experiment involves the 3d mixed-problem on 
an outer spherical shell, analogous to the outer annulus problem studied in 
Sec. 4.2.2. The second experiment is a crude two-domain model consisting of 
an inner cube (lying within the elliptic region and on which we place explicit 
sources) and an outer spherical shell (intersecting the type-changing cylinder). 
This single cube serves as a crude substitute for all the subdomains depicted 
in Fig. 2 (as well as the inner spherical shell not depicted). We would prefer 
an experiment in which a single cylinder replaced all these subdomains, but as 
yet have not developed a library of spectral routines for solid cylinders. Our 
simple two-domain model indicates that iteration between an inner elliptic 
region and an outer spherical shell is stable. Such an iteration also appears to 
be quite stable in the 2d setting. 
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4-3.1 Mixed problem on outer shell 



For the first experiment we define the outer spherical shell via A = e < r < 
R = 50, and consider the following exact series solution to the homogeneous 
3d HRWE: 



^p{r, e,ip) = Y^ ceoYeo{9, 0)-—- + XI XI '^irnYim{9, 0)ji{mna) 



1=0 



=1 m=l 



X ^j£{mQ,r) sm{m(p) + ni{mVtr) cos{mip) 



(91) 



here with r > a = 1 and Q = 0.2. The Qm are random expansion coefficients, 
obeying \cem\ < 1 and drawn from a uniform distribution, Y^rn(^, ^) is a scalar 
spherical harmonic, and je{z) and ni{z) are spherical Bessel functions. Notice 
that e < l^^l"^, whence the inner surface of the shell lies within the elliptic 
region, as discussed in Sec. 2.2. The form of this series is motivated by the 
rather more physical series 18 



2ky: 



-YUn/2,0)Y,o{9, 0)- 



oo 



,£+1 



2i+l 

+ 4KnY. E rnY,^{n/2,0)Y,^{9,0) 

1=2 m=2,4,6,--- 

X jiimVta) |^j£(mf2r) sin(m(y9) + niimVtr^ cos{mip) 

corresponding to two equal point charges, that is 

.6{r — a) 



6{cos6) 6{(p) + 6{(p + vr) 



(92) 



(93) 



in Eq. (10). The choice here of like signs for the point sources would make 
the problem similar to the gravitational problem, with radiation dominated 
by the quadrupole mode. (In the 2d case we chose equal and opposite signs to 
avoid the bothersome logarithmic 2d monopole.) We could of course base our 
test on a truncation of (92), but will instead work with (91) in order to test 
all modes defined for a given choice of -^max- 

We use the series (91) to seed inner boundary conditions at r = e, which 
would also be possible for (92), since the locations r = a, 6 = 7r/2, = 0, vr 
do not lie in the shell, and so the series is convergent for r > a. This requires 
a spherical-harmonic decomposition, which we perform with the NCAR rou- 
tines SPHEREPACK by Adams and Swarztrauber 47j . For this example, we 
have coded our spectral representation of the 3d HRWE as a matrix-vector 
multiply. As we do not explicitly form the matrix, we solve the resulting linear 



system using the iterative method GMRES[48|]. Our matrix-vector multiply 



involves integration preconditioning in the radial variable only. We use the se- 
ries above with i^^^ = 7, which corresponds to 1+3+5+7+9+11 + 13+15 = 64 
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modes. The capacity to perform spectral analysis and synthesis (transform and 
inverse transform) then requires that a corresponding physical grid has an an- 
gular resolution of at least Nq x A^^^ = 8 x 15 = 120 points. In anticipation of 
using such analysis and synthesis when dealing with nonlinearities, we there- 
fore work with a set of unknowns corresponding to Nj-xNgX = 60x8x15 = 
7200, although 60 x (120 — 64) = 3360 should be zero. Since only the remain- 
ing 60 X 64 = 3840 unknowns should be nonzero, as our linear solver we use 
GMRES, taking the number of iterations Mitr < 3840. As an example, with 
Mitr = 2000 we produce a numerical solution from which we plot and analyze 
various cross sections. Over the equatorial plane the numerical solution has a 
supremum error of 1.194 x 10~^. For Mitr = 1600 the corresponding supremum 
error is 1.038 x 10"^ 



4-3.2 Mixed problem on two coupled domains 



As the second experiment, we consider essentially the same problem described 
in Ref. llj, that is the linear inhomogeneous 3d HRWE of (10) with two 
Gaussian sources. 



g{r,e,Lp) = Qe 



(x-XH)'^+{y-yH)^+z^ 



(94) 



Our parameter choices are Q = 100, k = 3, with the location of the "hole" 
(source) at {xh.Vh) = (1,0) in the z = plane. Furthermore, we choose 
/i = —0.25 in order to ensure that the stability of our scheme is not predicated 
on symmetry. We determine the HRWE operator hj Q = 0.2. Our two basic 
sub domains are an inner cube and an outer shell, respectively determined by 

[^^minj^max] = [Z/min, l/max] = [^rnin; ^max] = 3] and [e, R] = [2.8,100]. 

To solve this problem numerically, we work with the two domains as decoupled. 
On the cube we use a pseudospectral method (point values as unknowns) based 
on GMRES, with the derivatives in the HRWE operator coded as a matrix- 
vector multiply. We currently employ no preconditioning whatsoever on the 
cube in our simple model, and so are limited in what resolutions we can obtain. 
For the experiment we have taken a A^a, x A^^ x A^^ = 30 x 30 x 30 truncation. 
For the outer shell, we again use GMRES but with the purely spectral method 
just described in Sec. 4.3.1. Therefore, we have adopted a hybrid pseudospec- 
tral/spectral approach. Our resolution in the annulus corresponds to a grid 
size of Nr X Ng X N^p = 62 x 5 x 9. Our GMRES-based experiment requires 
less than 5 percent of the combined storage required to perform an individual 
direct solve for each domain, although it does not achieve high accuracy. 

Starting with the cube, we iterate between the solves. That is to say, we first 
solve the problem on the cube with Dirichlet-zero boundary conditions. Via 
the global interpolation then available from the cube solution, we seed inner 
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Fig. 8. Cross sections of the 3d numerical solution. The top plot depicts 
the equatorial cross section on the outer shell, with multiplication by r enhancing 
1 /r fall off. The bottom plot depicts the equatorial cross section on the inner cube. 

Dirichlet data ^/^l^ for the shell solve. Since we are solving the 3d HRWE on 
a closed ball (cube plus shell) with no inner boundaries, the solution is only 
determined up to a free constant. We fix this constant by setting to zero the 
£ = mode of t/'Ie. The solution on the shell then provides corrected Dirichlet 
data for another cube solve, and so on. The iteration proceeds until the rms 
error between successive seed data sets ^/'le ceases to decrease (about 60 itera- 
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tions). Figure 8 shows the resulting numerical solution. For this low resolution 
experiment, the depicted outer shell solution has an absolute supremum error 
(as measured against a second numerical solution obtained with slightly more 
resolution on both the cube and shell, but with the same number of angular 
modes) measuring 1.478 x 10^^ over the equatorial plane. Although we have 
not considered large values of Q, we note that the iteration appears to be com- 
pletely stable and convergent. An instability could set in at higher resolutions. 
Nevertheless, we find the results from this model 3d code as well as similar 
iteration in the 2d setting encouraging. 



5 Summary, conclusions, and outlook 

The motivation for this work arose in problems of gravitational waves and 
black holes. In a late stage of inspiral the interaction of the holes can be 
strongly affected by nonlinearities, while the radiation reaction due to gravi- 
tational wave emission is weak. It is useful to approximate such a system as 
having helical symmetry, in which the sources and fields rotate rigidly. Two 
aspects of this physical problem are then important to the mathematical meth- 
ods explored in this paper: (i) "Helical reduction," that is, imposing helical 
symmetry, converts the PDEs describing the fields into a system of helically 
reduced wave equations (HRWEs), a mixed system that is elliptic inside, and 
hyperbolic outside a type changing surface, (ii) Nonlinearities are significant 
only interior to this type changing surface. In this paper we have divided the 
topologically nontrivial domain of the problem into basic subdomains. Cru- 
cially, only a single subdomain contains the type changing surface and the 
hyperbolic region exterior to it, the outer annulus/shell subdomain (annulus 
for 2d, shell for 3d). The other subdomain or subdomains span an elliptic 
problem. We have therefore confined aspect (i) to a topologically simple outer 
subdomain, and aspect (ii) to the subdomains spanning an elliptic problem. 

Sections 2 and 3 of this paper have focused on the model problem of the 
2d HRWE for a linear scalar field, though some generalization to 3d and to 
nonlinear problems has been included. Likewise, while the numerical results 
in Sec. 4 have predominately focused on the 2d linear scalar problem, we have 
also included some results for a nonlinear 2d scalar model and for a linear 3d 
scalar model. To obtain our numerical results we have necessarily made many 
choices in the details. This was also case in Sec. 3; clarity of presentation 
required some specificity in the adopted methods. It is, however, important 
to distinguish the full set of choices made from what we believe are robust 
conclusions that transcend many of the choices. In the present section, we 
therefore remark on what, in our view, are the central lessons of this work. 

The choice to use a spectral (or a pseudospectral) method does not need 
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justification for an elliptic problem. Our mixed problem inherits much of that 
justification. In particular, we are solving a boundary value problem, so a large 
set of equations must be solved simultaneously. For such problems spectral 
methods, which reduce the number of unknowns, are particularly valuable 
and widely used. A less general issue is the nature of the decomposition of 
the two center domain of the problem into subdomains. The choice illustrated 
in Fig. 1 is very much driven by the convenience of rectangular and annular 
subdomains, for on such regions spectral methods are well developed. For 
analogous reasons, the similar decomposition in Fig. 2 would be used for the 
3d problem. 

What we wish to emphasize here is that one should make a distinction between 
the choices for the outer annulus/shell and those for the inner region, the sub- 
domain or concatenation of subdomains for the elliptical region. (One could 
also divide the outer shell/annulus into further concentric sub-shells/sub- 
annuli, but this possibility only adds further technical complexity). It is quite 
conceivable that the optimal choice for the outer region is not the same as that 
for the inner region. One might, for example, treat the inner elliptic problem 
by whatever method most conveniently or efficiently generates a set of equa- 
tions for the interior, and then use the spectral method of this paper only for 
the outer region with its mixed PDEs. Even more freedom is available if the 
inner and outer regions are solved separately, with an iteration between them, 
as was done in the second numerical experiment of Sec. 4.3. In this case, one 
could use separate iterations to solve within each region, each within the outer 
loop iteration that feeds values from the inner to the outer region and vice 
versa. We will return to this possibility in a moment. 



We have made the choice of using integration preconditioning (IPC) in all 
subdomains. Whether or not this is an optimal choice remains unresolved. IPC 
is known to be efficient in reducing the condition number of the coefficient 
matrix for certain types of ODE boundary value problems. Our study of a 
novel boundary value problem, one heretofore not considered in the context of 
spectral methods, has to some extent further elucidated the utility of IPC in 
such ODE settings. The results in Table 4 demonstrate that IPC reduces the 
condition number of the matrix stemming from a problem not considered in 



our primary reference [40[]. Moreover, our work demonstrates that IPC affords 
a relatively direct way to handle a multidomain scenario in two dimensions. 
(See also 45|] for an application of IPC in a Id multidomain setting.) However, 
it remains to be seen if its advantages will carry over to 3d and to PDEs 
significantly different from those of our test problems. 



Coming to grips with whether or not IPC is an appropriate method for 3d PSW 
problems, requires that we ask some broader background questions. In partic- 
ular, is it important to reduce the condition number of the coefficient matrix? 
In Sec. 4, solutions to 2d problems were found by straightforward Gaussian 
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elimination for which conditioning issues are not as critical. For 3d problems, 
however, we almost certainly expect to use an iterative Krylov method, and 
GMRES in particular. The efficiency of such iterative solvers is known to de- 
pend sensitively on condition number. But such methods are also sensitive 
to other aspects of the coefficient matrix, in particular to the distribution 



of eigenvalues [461]. In any case, owing to the relevant domain decomposition 
(Fig. 2) for a 3d problem, insofar as the elliptic problem is concerned the rel- 
evant issue is the effectiveness of IPC for the 3d HRWE on cylindrical and 
(inner) spherical-shell domains. As noted in [i^, S^], special alternate tech- 
niques are often necessary for polar and cylindrical subdomains which contain 
a center or central axis. 



We turn now from the open question of whether IPC could be effectively 
applied to all, or almost all, subdomains of a 3d PSW problem, to the rather 
clearer question of IPC for the outer region itself. For a scheme in which inner 
and outer regions are separately solved, we have some enthusiasm for the 
advantages of IPC. Our enthusiasm stems from the fact that in the outer region 
the Fourier, or spherical harmonic modes do not mix for a linear problem. 
This nonmixing originates in the symmetry of the problem, the fact that the 
annulus or shell is concentric with the coordinate surfaces of constant radius. 
The situation is closely analogous to the reason that separation of variables 
applied to this region leads to a set of decoupled radial ODEs. In line with 



theoretical and numerical arguments presented in [40|, our studies indicate 
that IPC is a promising way of solving the system of ODEs arising in the 
outer region of our problem. 



These simple considerations do not apply, of course, if the problem is nonlinear. 
The nonlinearity will mix Fourier or spherical harmonic modes. However, the 
nonlinearities in the underlying physical problem are very weak in the outer 
region. (Here we are assuming that the inner surface of the outer region is 
chosen to be near the maximum radius that allows it to be inside the type 
changing surface and to have the needed overlap with the subdomains of the 
interior region.) As a result, there would seem to be number of possibilities 
for handling the nonlinearity. First, it is quite plausible that it is adequate to 
treat the problem as linear in the outer region. If it is not adequate to treat 
the nonlinearity as null, then one could exploit the fact that it is weak. The 
weakly nonlinear problem could be solved by iteration in which each iteration 
generates a "known" source term, and that source term is projected onto the 
angular mode basis before the system is solved. Just as would be the case 
with separation of variables, this procedure would not mix modes (after the 
projection of the known source is done) so the solution would again be highly 
efficient. The weakness of the source should guarantee that a few iterations 
would be sufficient to get the needed resolution. Finally, a Newton-Krylov 
method might be based upon use of the outer domain's spectral transform and 
inverse transform (with all nonlinear operations performed in point space). For 
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a weak nonlinearity the associated Jacobian should not significantly perturb 
the basic linear operator. 



Our next steps with these methods will parallel those of Refs. |8|- 19|. The 
method used in those studies is, by its nature, incapable of high precision, 
whereas our approach is capable of the precision that is needed for some 
purposes (initial data for evolution codes and studies of radiation reaction for 
instance). In particular, we will next formulate 3d linear and nonlinear scalar 
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model HRWE problems. We will then use the techniques of Ref. 
linearized gravity, and those of Ref. 19|] to do gravity with a post-Minkowski 
approximation. Finally, we will treat the problem in fully nonlinear general 
relativity. 
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